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INTRODUCTION 



I . 

Ocean Thermal Energy Conversion (OTEC) plants will be 
either moored or dynamically positioned in the ocean and, 
accordingly, will be subject to wind, wave and current load- 
ing. This report deals with the motion of a floating body 
of arbitrary shape subject to wave motion. The interaction 
of the waves with the floating body is based on the inviscid 
(potential flow) theory. This assumption appears to be 
valid as long as the body involved is large compared to the 
amplitude of the relative fluid motion. 

This report describes the first phase of a project to 
analytically determine the dynamic response of large float- 
ing OTEC plants to ocean waves. The long range goal of the 
project is to develop analysis and computer code for compu- 
tation of the dynamic response and structural loading of 
OTEC structures of rather general configuration. Such 
configurations are considered to be composed of a large 
displacement parts where diffraction theory will be required 
and/or small appendages or members which can be dealt with 
by use of a Morison equation type analysis. This phase of 
the work deals with the analysis of large displacement 
bodies of arbitrary shape. The analysis is based on the 
use of a surface distribution of sources using quadr ala t er al 
source patches or panels. The source strengths over each 
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panel is assumed constant in this phase of the project. The 
next phase of the work will investigate the use of triangular 
panels where the source strength varies linearly between the 
values at the corner node points. A comparison of the tri- 
angular patch method will be made with the quadr ala t er al 
source patch method with constant source strengths over the 
patch. The better of the two methods will be used in the 
future phases of this work. 

In particular, this report deals with the analytical 
evaluation of the wave forces and overturning moment acting 
on a large bodyof arbitrary shapewithout appendages in water 
of finite depth. The same structure is also considered to 
oscillate in surge, heave, sway, roll, yaw and pitch and for 
such motion the added mass and damping coefficients are deter- 
mined. Then, the equations of motion for a free-floating 
body are developed, and using the computed wave excitation 
forces and moments and added mass and damping coefficients 
the response of the floating structure is evaluated. 

Although the problems of wave interaction with a fixed 
body and the oscillation of the same body in otherwise 
still fluid appear to be physically distinct, they are 
mathematically similar and as a result are dealt with 
simultaneously herein. The primary mathematical difference 
in these problems is the kinematic boundary condition which 
is applied on the immersed surface. In all cases, however, 
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the velocity of the fluid in the normal direction relative 
to the immersed surface must be zero but this statement takes 
a different form in each case. 

The method used to describe the fluid motion involves 
the use of a Green's function or source distribution. The 
rigid, immersed surface of the structure is represented by a 
surface distribution of sources which is assumed to be con- 
stant over quadrala teral panels and the zero relative normal 
velocity condition is applied in order to determine the 
strengths of the sources. Once the source strength distri- 
bution is, known the velocity potential at some point in the 
fluid region is determined by summing the effect of all of 
the sources at the point in question. 

Several authors have computed hydrodynamic coefficients 
associated with specific shapes. The added mass and damping 
coefficients for a s emi- immer s ed floating sphere undergoing 
heaving motion has been obtained by Havelock (2). The vertical 
wave excitation force was not calculated by Havelock but 
using the Haskind's Relations [Ref. 3] this quantity may be 
determined from the damping coefficient. Kim (4) determined 
the added-mass and damping coefficients for this same config- 
uration to wave excitation in heave and surge. Although 
calculated for the infinite depth case, these results pro- 
vide an important comparison for establishing confidence in 
the mathematical validity of the response calculations based 
on the present analysis. 
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Several authors have also considered the wave interaction 



with large displacement fixed and floating bodies. In the 
case of large North Sea Gravity platform [Refs. 10 and 20] 
the agreement between the theory and experiment is almost exact. 
In such structures the caisson was we 1 1- s ubme r ge d and the 
superstructure only pierced the free-surface . Thus, nonlinear 
effects which tend to be most pronounced near the free surface 
were absent. Additional experimental data and comparisons 
with the linear theory for the case of fixed bodies which 
are considered to be significant include those presented by 
Garrison and Seetharama Rao (6) for the case of a bottom 
mounted hemisphere. Unfortunately, the wave amplitudes were 
rather small. Also, Garrison and Chow (7) have presented a 
comparison of the theoretical results based on the distributed 
source theory similar to that developed herein with two slightly 
different submerged oil storage tank models. Unfortunately, 
the vertical force part of this data was rather poor because 
the vertical model supports had very low spring rates making 
the natural frequency of the model support system too close to 
the wave excitation frequency. It is significant, however, 
that these tests represented North Sea design conditions and 
in the case of the horizontal force quite good agreement was 
obtained between the linear theory and experimental results. 

A third experimental study involving a short vertical circular 
cylinder (1.16 cylinder radii water depth) extending from the 
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bottom through the free surface should be noted. Chakrabarti (8) 
has made measurements of the horizontal force and overturning 
moment for this configuration and found the results to compare 
well with MacCamy and Fuch T s (1) corresponding theoretical re- 
sults. However, Chakrabarti gives no indication of the wave 
height involved in Ref. (8). It is noted that the MacCamy and 
Fuchs solution is simply a special case of the present more 
general method so that agreement with this solution implies 
agreement with the present method. 

A series of tests have been conducted by Hogben and 
Standing (9) of the National Physical Laboratory, London, 

England with fixed vertical cylinders of square and circular 
cross-section. They also made calculations using a computer 
program based on diffraction theory similar to the one pre- 
sented here and good agreement between theory and experiment 
has been observed. 

Only limited experimental data has been presented in the 
literature for the dynamic response of three-dimensional 
floating structures to waves. Among these the results of 
Faltinsen and Michelsen (11) are the most complete. They 
compared calculations of diffraction theory with the heave 
response of a floating box and the agreement obtained was 
quite good. Also Garrison (12) showed excellent agreement 
between the results of diffraction theory and the heave and 
pitch response of a disc buoy. However, as a general assess- 
ment of the experimental data available in the literature, 
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for the hydrodynamic coefficients and dynamic response of 
floating bodies, it appears that the amount is very small and 
the understanding of viscous effects and differences or sim- 
ilarities of the experimental results with the theory is non- 
existent. This is an area where additional work is needed. 

2 . FORMULATION OF PROBLEM 

Consider a rigid object of arbitrary shape and having 
characteristic lineal dimension a with center of gravity 
submerged to a depth d beneath the free surface in water of 
depth h as shown in Fig. 1. The structure is considered to 
be smooth to the extent that its unit normal vector is a 
continuous function and it may or may not intersect the bottom 
or free surface. Two coordinate systems are identified, x,y,z 

coordinates with origin fixed at the free surface and the body 

? ? ? _ _ 

coordinates x,y,z coordinates located at depth y = -d . The 

bars over the symbols denote dimensional quantities. 

The mathematical problem which is now established is that 
associated with the fluid motion, pressures and resulting forces 
induced by the small amplitude oscillation of the object in 
its six degrees of freedom as well as the fluid motion asso- 
ciated with the interaction of the fixed object with a train 
of regular waves. The small amplitude oscillatory motion of 
the structure about its equilibrium position with frequency a 
is described by the relationship, 

X. (t) = X? R [e‘ i0t ], i = 1,2,3 (2.1a) 

i i e 
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FIGURE 1 



DEFINITIONS 



0 i (t) = 0° R e te 10t ], i = 4,5,6 



(2.1b) 



where the subscript i = 1,2,3 denotes lineal oscillation in 

_!_!_! 

the x,y,z directions, respectively, and i = 4,5,6 denote angular 

I I I 

oscillation about the x, y,z axes, respectively. The real 



numbers X° and 0° denote the amplitudes of the motion. The 

displacements, X^, and X^ are normally referred to as surge, 

heave and sway, respectively, while the angular components 0^ , 0,. 

and 0^ are called roll, yaw and pitch, respectively. The second 

problem dealt with simultaneously is the interaction of a train 

of regular surface waves with the object fixed in space. The 

regular incident waves of wave height 2^° = H and wave length 

* 

L are assumed to progress in the positive x-direction and 
interact with the fixed structure. 

The fluid is assumed to be incompressible, and the motion 
irrotational and harmonic with time dependence e iat in all 
cases. It follows, therefore, that a velocity potential exists 
such that the fluid velocity vector may be defined as 



q- = R £ fV(|> j (x,y , z) e 1CTt ], j = 1,2, ..6 



( 2 . 2 ) 



where [ 4> ^ e iat ] denotes the velocity potential associated 

with the motion induced by oscillations in the six degrees of 



freedom and V = i — — + j — 



+ k , the symbols i,j,k denoting 



8 x 8y 8 z 

the unit vectors in the x,y,z directions, respectively. 

For the case of regular wave interaction with the fixed 
structure, the velocity potential may be written as the sum 

1 = o + <t> 7 



(2.3) 
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where <f> 0 denotes the potential associated with the incident 
wave in the absence of the structure and <j> ^ denotes the 
scattering potential which results from the presence of the 
structure in the wave train. In this case, the fluid velocity 
vector is given by 

-* -> 

q' = R e [ V ( 4 , 0 + <f, 7 )e -iat ] (2.4) 

The continuity equation shows that <(>_., (j = 0,1,2, ...7) 
must satisfy the Laplace equation 

v 2 4 > (x, y , z) = 0 (2.5) 

and from the linearized form of Bernoulli’s equation, which 
is applied throughout, the dynamic fluid pressure is given by, 



Pj = R e [i p a e i<Jt ], j = 1,2, ...6 (2.6) 

For the second problem involving wave interaction with 
the fixed object, the pressure is given by 



P' = R [i pa ( 4 o + <}>-,) e 10t ] (2.7) 

e / 

The velocity potentials must satisfy certain boundary con- 
ditions in addition to Eq. (2.5). These include the linearized 
free surface boundary condition. 



— — j ( x , 0 , z ) ■ § $ j (X , 0, z) = 0 , j = 0,1,2, ...7 

3y 



( 2 . 8 ) 
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Equation (2.8) is familiar from linear wave theory. Also, 
must satisfy the kinematic boundary condition on the bottom in 
all cases , 

3j> (x,-h,i) =0, j = 0,1,2 7 (2.9) 

8y 



The surface of the structure in its mean position is 
described by 

S( x , y , z) = 0 (2.10) 

When the structure oscillates, the velocity of the fluid in 
the direction normal to the immersed surface must equal the 
velocity of the surface normal to itself. Also, for the case 
of wave interaction with the fixed body it is necessary that 
the normal velocity component be zero. These conditions are 
stated mathematically as 



9 h 

a ? 

9 h 

9_^_3 

9 n 

4 

9n 

115 

9 n 

1^6 

9 n 

d± 7 

9 n 



o 



- ia X. n 

1 X 




(a) 


O 






-ia X 0 n 

2 y 




(b) 


_ o 






-ioX.n 
3 z 




(c) 


O 






-iaG^ [ (d+y) 


n - zn 
z y 


(d) 


O 






-io0, [zn 
5 x 


xn 

z 


(e) 


o 

- i a0 , [ xn 
6 y 


(d+y) n x ] 


(f) 


O 

^lie 

ro | 

1 




(g) 



( 2 . 11 ) 
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where n = in + in + kn denotes the unit normal vector on 
x y z 

the surface of the object directed outward into the fluid and 
d denotes the depth of submergence of the x, f y ’z ’coordinate 
origin. Finally, the velocity potentials must satisfy the 
radiation condition which allows only outgoing waves. 



.(? 1 6,;;)->.(9)^-‘ 1 cosh[k(Y + hlie lkl i - j. . 1,2, ..7 

cosh(kh) 1 



(2 



-2 -2 h -1 - - 

where r^ = [x + z ] and 0 = tan (z/x) . The wave number is 

defined as k = 2 tt/L where L denotes the wave length and 

is related to the frequency of the motion according to the 

well-known relationship, 

2 



— = k tanh (kh) 

g 



(2.13) 



The velocity potential of the incident wave alone progres- 
sing in the positive x direction which satisfies Eqs . (2.5), 

(2.8) and (2.9) is given by 



_ o 



<f> 0 ( x , y ) = 



i g n .,. cosh[k(h+y) ]e 
a cosh (kK) 



ik(x cos ^ + z sin $) 



(2.14) 



where as indicated in Figure 1, ^ denotes the angle of the 

- o — . 

incident wave measured clockwise from the x-axis and f|* = H/2 
denotes the amplitude of the incident wave, H being the wave 
height. Moreover, from the linearized form of Bernoulli’s 
equation, the free surface of the incident wave alone is 
evaluated from 

1 ^ ( — i n t \ 

(2.15) 






. 12 ) 
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Using (2.14) this gives, 

— ° 

% = ^ cos(kx cos f) + kz sin i|j) (2.16) 



Thus, for purposes of reckoning phase angle it is noted that 
at t = 0 the crest of the incident wave is just passing the 
coordinate origin. 

For convenience in carrying out the solution for the seven 
potentials, <j> , and to show clearly the dependence of the 
solution on the parameter, a = 27ra/L = ka, the relative water 
depth, h = h/a and the relative depth of submergence, d = d/a, 
the space variables and amplitudes are first made dimensionless 
with the characteristic linear dimension of the object, a. 



x = x/a, y = y/a, z = z/a, r = r/a 

o o _ o o 

X i = X i /a ’ (i = i’ 2 * 3 )’ (i = 4,5,6) 



° _° _ _ 2 — 

n* = n /a = H/2a, r ^ = v^/a, V = a a/g 



(2.17) 



and then the dimensionless potential functions u^ are introduced 
as , 



io <f>j ( x , y , z) / gaX_. = a tanh(ah) u_.(x,y,z), j = 1,2, ..6 

° 

ic <f> 7 (x, y , z) /ga n* = -a u 7 (x,y,z) 

O 

io <f> 0 (x, y) /ga = -a u Q (x,y) 



( 2 .18a) 
( 2 . 18b) 
(2.18c) 



The complex dynamic pressure amplitude can now be written 
by use of the linearized form of Bernouli's equation (2.6) and 
(2.7), as 



19 



Pj = a tanh (ah) u (x,y,z), j = 1,2, ...6 



(2.19a) 



P' 



, r /-uj_ \i ia(x cosHz sinip) , . .. . 

_ cosh [a (h+y ) ] e r -a u^(x,y,z) (2.19b) 

cosh (ah) 



where the complex amplitudes of the pressure, p^, ( are defined 
P . 



as 



= [p . (x,y , z) e 1CTt ], j = 1,2, ...6 



p gaX“ e J 



(2.20a) 






(2 . 20b) 



The boundary value problem which describes the fluid motion 
arising from the oscillation of the rigid object in its six 
degrees of freedom as well as the scattering of the incident 
wave may now be written concisely in terms of dimensionless 
parameters. The potentials u^(x,y,z), i = 1,2, ...7, continuous 
in the fluid region is sought such that: 



V Uj (x,y , z) = 0 

9_u . (x , o , z) -a tanh(ah) u.(x,o,z) = 0 

By J J 

3 u . ( x , -h , z) = 0 

By 3 



9u_.(x,y,z) = g . ( x , y , z) on S(x,y,z) = 0 

9n J 2 



, A x „ x -h cosh [ a (h+y ) ] e iar 1 

U (r,0,y)-X 1 (e)r 1 ^thUh) *°> r l 



( 2 . 21a) 
(2 . 21b) 

(2.21c) 
( 2 . 21d) 

( 2 . 21e) 
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The g ( x , y , z ) functions represent the prescribed functions which 
depend on the mode of oscillation (j = 1 , 2 , ... 6 ), and j = 7 
corresponds to scattering. These functions represent the dimen- 
sionless form of ( 2 . 11 ) and are given by 

g l = n x’ g 2 = n y’ g 3 = n z ( 2 . 2 1 f ) 

g 4 = (d+y)n z -zn y , g 5 = zn^xn^ g & = xn y -(d+y)n x 

g 7 = cosh (ah) ^ n y s 1 nh [ a (h+Y > 1 + i cosh[a(h+y)] 

/ llt , . in -i i a ( x cos ^ - sin ^ ) 

.(ncosV+nsinV9Je 
x z 

Equations (2.21) define seven b ounda r y - va lue problems which 
correspond to oscillation of the immersed surface in each of 
the six degrees of freedom and to scattering of the incident 
wave by the fixed body. The problem statement for each of the 
potentials is the same; the only difference lies in the form of 
the boundary condition to be applied on the immersed surface. 

3. REPRESENTATION OF THE POTENTIAL 

The boundary value problem for oscillation in the six degrees 
of freedom and scattering of the incident wave is specified in 
(2.21). The solution, (i.e.), the function u ( x ,y,z), may be 
represented by use of a Green's function having the physical 
interpretation of a point wave source of unit strength. These 
sources are distributed over the surface of the object according 
to the source strength function, f, so that the potential at 
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some point (x,y,z) within the fluid region is given by the 
surface integral 



Uj(x,y,z) = fj(5»n,C) G (x , y , z ; C , 0 , O dS (3.1) 

where ( £ , n > O represents a point on the surface of the struc- 
ture, G denotes the Green's function (or source potential) and 

- _ 2 

dS = dS/a denotes the dimensionless surface area element. The 
Green's function is defined, therefore, as the function which 
satisfies 

V 2 G (x, y , z ; £ , n , C) = S(x-£) 6(y~n) S(z-0 (3.2) 



as well as the boundary conditions (2.21b, c, and e) . Such 
a function is given by Wehausen and Laitone (13) as 

G(x,y , z ; £ , n , 5) = ^ + G* (3.3a) 

where G* = ^ 

R 



+ 



2P.V. 




( y+v) e yh cosh[y (n+h) ] cosh[y(y+h)] 
y sinh (yh) -v cosh (yh) 



J 0 (y r) dy 



+ i 



2 2 

2 tt (a -v ) 



cosh [a ( n+h) ] 

2 2 
a h - v h 4- v 



cosh 



[a (y+h) ] J n (ar) 



(3.3b) 



R = [(x-O 2 


+ 


(y-n) 2 + ( z- <;) 2 ] 


(3.3c) 


R’= [(x-5) 2 


+ 


( y+2h+n ) 2 + (z-c) 2 ] -5 


(3.3d) 


r = [ (x-C) 2 


+ 


(z-O 2 ] 1 " 


( 3 . 3e) 


2- , 

v = a a/ g = 


a 


tanh ( ah ) 


(3 . 3f ) 
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and P.V. indicates principal value of the integral. An alter- 
nate series form of the Green's function is also given by 
Wehausen and Laitone (13) as 

2 2 

2 7T ( V ) 

G ( x , y , z ; £ , n , <; ) = z cosh [ a (h+y ) ] co sh [ a (h + n ) ] [ Y (ar)-iJ (ar) ] 

a h-v h+v ° ° 



+ 4 



2 



o-k ♦ v 2 > 



(3.4a) 



k=l M k h+ V h_V 



cos[u k (y+h) ]cos[n k (n+h) ] K q ( u k r) 



where J q and denote, respectively, Bessel functions of the 

first and second kind of order zero and K denotes the modified 

o 

Bessel function of the second kind of order zero. The quanti- 
ties are the real positive roots of the equation 

u k tan (p k h ) + v = 0 (3.4b) 

The Green's functions specified in (3.3) and (3.4) are 
equivalent. However, for purposes of numerical evaluation it 
is found that the series formulation given in (3.4) can be 
evaluated more accurately and with much less expenditure of 
computer time than the integral form when (ar) is not too small. 
On the other hand, for small and zero values of (ar) the integral 
form given in (3.3) is necessary. The integral form of the 
Green’s function as given in Eq . (3.3) may be integrated 

directly or, it may be rearranged as in Appendix C so as to be 
somewhat better conditioned for numerical evaluation. In the 
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computer program when the integral form of the Green’s function 
is required, the form given in Appendix C is actually used for 



numerical evaluation. It appears that it is somewhat less time 
consuming to evaluate than the form given by Eq . (3.3) . The 

potentials u^ are represented by (3.1) and with the introduc- 
tion of the Green’s functions (3.3) and (3.4) (or equivalent form 
given in Appendix C) the only unknown in (3.1) is the source 
strength function f ^ . To evaluate this it is necessary to 
take the derivative of (3.1) in the direction normal to the 
immersed surface and then apply the boundary condition (2.21d). 
This results in the following integral equation from which f ^ 
is to be determined: 



I f/” f (r n r N iG(x,y , z ; c , n , C) ds = g.(x,y 



, z), j = 1,2. ..7 



(3, 



where 3 G/ 3n is obtained using the form VG . n where VG is deter- 
mined by straightforward differentiation of (3.3) or (3.4). 

The normal derivative of (3.3) is 



3 G 

3 n 



= - — 3 [n (x-O + n (y-n) + n (z-c) ] 



— , 3 [n (x-£) + n (y + 2h + n) + n ( z - C ) ] 

iv X V z 



- 2 P . V 



j ~ » p(p+v)e ^fcoshlu (y+h)1 



M sinh(yh) - v cosh(yh) 



(3.6) 



cosh [ n (y+h) ] J. (pr) 



[n (x-O+n (z-?)j-n s inh [ p (y+h) ] J ( yir ) 
x z y ° J 



27Ta(a 2 -v 2 )cosh[a(ri + h) ] T cosh[a(h+y) ] J^(ar) 

-i j 2 L 

a h - v h + v r 

-n sinh [ a ( h+y ) ] J ^ ( a r ) 

and the normal derivative of the alternate rearranged form of 
Eq . (3.3) is given in Appendix C. The normal derivative of 

the series form given in (3.4) is 



(n x (x-C)+n z (z-c) ) 



2 2 

3G 2 tt ( v -a ) a cosh [ a (h + n ) ] T . , r ,, . s , r , x , x i 

- 2 2 s i nh [ a ( h+y ) ][f (ar)-iJ (ar)jn^ 

ah-vh + v L 



“ [ n x ( x "5) + n z ( z “ O] cosh [ a (h+y ) ]~ [y ± ( ar ) - i ( ar )] 



(3.7) 



°o o 2 

E u k (y^ +V ) cos [M k (n+h) ] p 

.. 2 „ „ .. 2 U — : [ si 



k-i u k h + v2h - v 



sin[ y k (y+h) ] K c ( M k r) n y 



K 1 ^ u b- r ) 

+ cos[u k (y+h)] [n x (x-£) + n ? (z~c)] 



where n , n and n denote the components of the outward unit 
x y z 

normal vector on the immersed surface. 

Returning to the integral equation (3.5), it is evident 

that the 3G/3n occurring therein as given by (3.6) is singular 
3 

like 1/R as the point ( £ , n , £ ) approaches the point (x,y,z). 
Thus, special care must be taken when integrating this term 
of (3.6) over the singularity. Moreover, the kernel of (3.1) 
is singular like 1/R at this same point and likewise special 
attention must be given to this situation. Taking the integral 
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over the complete surface of the object as the sum of two parts, 
one part being a small circle of radius r^ about the point 
(x,y,z), as indicated in Figure 2, and the second part being 
the remainder, S ? , of the surface, we may write (3.1) as 



hffi V c>n> ° If <x ' y ' z;5 ' n ’ ,:) dS 



+ 




f j U > n > O 



8 G , 

(x * y * z; 



K , n , O dS 



gj (x ,y , z) 



(3.8) 



where £ denotes the area inside the small circle of radius r . 

o 

Now, using the expression for 3G/8n as obtained from (3.3a) 
and keeping in mind that f (£,n>C) is a well-behaved function, 
we may take the limiting case as £, or equivalently, r Q tends 
to zero to obtain 



1 im 
I->0 



1 

4 TT 



f (x , y , z) 





dS 



+ 



1 im 
l-y 0 



47 £ j (x,y ' z) 




8 G * 

9 n 



dS 



(3.9) 



dim 1 /7* , , . x (x,y,z;£,n,0 dS = g.(x,y,z) 

Z-+0 4 TT //s' r j u,n ’ 9n J 



Letting the point (x,y,z) lie along the normal at a distance 
e above the surface as indicated in Figure 2, it is evident 
that the first integral in (3.8) can be expressed as 
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X! 






_C 




2 7 



FIG. 2 SINGULAR POINT 



1 im 

e+0 



h f j (x ’ y ’ z) ./^k ( R )dS = lim n h f j (x ’ y ’ z) lr/o 0 



r ->0 
o 

e-> 0 



1 im 1 

= r o' >0 [ 47 V x,y,z) 2n ( 7=f=; 

e +0 /e +r o 



2ttt dr 

/ r ^ + e z 

(3.10) 



•1) ] = - - f j (x,y , z) 



provided is chosen to be small enough that the surface area 

inside r^ may be considered to be plane. Moreover, since G* 
occurring in (3.9) is regular, the second integral vanishes 
in the limit as giving the integral equation 



-f.(x,y,z) + A ff (X ’ y ’ Z! - 2 ® j (x,y,z) (3.11) 



In (3.11) the singular point at (x,y,z) is considered to be 
excluded from the surface S. 

As noted, in (3.1) the Green’s function is also singular 
and special consideration is needed. However, using a procedure 
similar to the above, it can be shown that this 1/R singularity 
contributes nothing to the surface integral in (3.1). Thus, 
(3.1) may stand as is without modification and the point (x,y,z) 
considered to be excluded from the surface S. 

However, the integral of 1/R over the facet of area AS 
outside the singular point is not zero. This will be discussed 
in detail in Appendix B. 
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4. 



HYDRODYNAMIC FORCES AND MOMENTS 



The forces and moments caused by the dynamic fluid pressure 
acting upon the immersed surface of the structure may be ob- 
tained from the integrals. 



F i j(t) = -(1.0 or a)^/^Pjg i dS, i,j = 1,2, ...6 
F^(t) = -(1.0 or a )J‘J§ P ' g ^d S , i = 1,2, ...6 



(4.1a) 



(4.1b) 



where F (t) denotes the i-th component of load arising from 
the j-th component of motion and F^(t) denotes the i-th compon- 
ent of wave force (or moment). The coefficient 1.0 is to be 
used in the case of a force (i = 1,2,3) while a is to be used 
when F denotes a moment (i = 4,5,6). The sign convention of 
the forces, moments, displacements, velocities, etc. follow 
the right hand rule. 

It is convenient and conventional to place the forces and 
moments in dimensionless coefficient form. Accordingly, the 
following complex wave force (and moment) coefficients are 
defined using (4.1b): 

F . ^ i 

C. = — ■■ , i = 1,2,3 (force) (4.2a) 

1 O o 

pga n 

* 

F . . , e 16 i 

C. = — — ^ - ax; , i = 4,5,6 (moment) (4.2b) 

pga n 

* 
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The symbol F_. (max) denotes the maximum value of the oscillatory 
force (or moment) and is taken as positive. The wave force 
may, therefore, be defined by the complex number or by the 

magnitude of and the phase angle 6_^. Once | C _^ | and 6^ are 

known, the wave force (or moment) is expressed as a function 
of time as 

F^(t) = (1.0 or a)p ga^ (H/ 2a) | | cos (6^ - at) (4.3) 

where the coefficient (1.0) is applicable in the case of a 
force (i = 1,2,3) and (a) is applied when (i = 4,5,6) and F (t) 
denotes a moment. The dimensionless wave amplitude is defined 

o 

as n* = H/2a. It is further recalled that all phase angles 
are measured in relation to the incident wave; at t = 0 the 
crest of the undisturbed incident wave is located at the coor- 
dina te origin. 

According to the definitions (4.2), may be expressed, 

using (4.1b) with (2.20b) and (2.19a), as 



C . 

l 




u 7 ( x , y , z) 



cosh [a ( h + y ) ] ia(x cos i/j +z sin^) 
cosh (ah) e 

(4.4) 



.g jL (x,y,z) dS, i = 1,2, ...6 

Thus, once u^(x,y,z) at all points on the immersed surface is 
determined, the complex wave force (or moment) coefficient 
may be evaluated by evaluating the surface integral indicated 
in (4.4). 
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Turning now to the force (or moment) produced by oscilla- 
tions in the various degrees of freedom, (2.19a) and (2.20a) 
may be used in conjunction with (4.1a) to obtain 



F i j ( t) =- ( 1 . 0 or a) pga 3 X° vR g [e 10 t ff s u j ( x , y , z ) g ± ( x , y , z ) d S ] , 

i , .1 = 1 , 2 . . 6 



(4.5) 



where as previously noted (1.0) corresponds to the case 
force (i = 1,2,3) and (a) corresponds to the moment (i 
Furthermore, if the two dimensionless real numbers M.. 
are defined as 



of a 



= 4,5,6) . 

and N . . 
ij 



M ij = -^//^^ j (x, y , Z ) g i (x, y ,z) dS 
N ij = -I">j^- j (x,y,z)g i (x,y,z) dS 



(4.6) 

(4.7) 



where Re and Im denote real part and imaginary part, respectively, 
then (4.5) may be written as, 



F..(t)=(1.0 or a) pga^X.vR [(M.,+iN..)e ^ at ], i,j=l,2,..6 (4.8) 

il 1 e ij ii 

where, as in (4.5), the coefficient (1.0) corresponds to i=l,2,3 
and (a) corresponds to i=4,5,6. 

Equation (4.8) may be further rearranged by defining the 
dimensionless parameters, 



F, . (t) 



C 1.U) 






, i = 1,2,3; j = 1,2, ...6 



P ga' 



and 



F. . (t) 



C ij (t) 



. ±3 



_4 
p ga 



, i = 4,5,6; j = 1,2, . . .6 



(4.9a) 



(4.9b) 
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2 - 

With (4.9) and using the definition v =a a/g,(4.8) may be 
written. 



C. .(t) 
ij 



--M. .X. (t)- — N. .X. (t) , i, j 
g ij J g 1J J 



1 , 2 , . . . 6 



(4.10) 



and N (i,j = 1,2, ...6) represent the inertia and 
damping tensors, respectively, having 36 elements each. The 
first index denotes the direction of the force or moment and 
the second index is associated with the component of the motion. 
It can be shown, moreover, that these tensors are symmetrical 
so that 



M. . = M. . 
iJ J i 



(4 .11a) 



and 



N . . = N . . 

J 1 



(4 . lib) 



Except for the diagonal terms it is rather difficult to 
give simple definitions for the dimensionless added mass and 
damping tensors, M _ and N _ , respectively. They are best 
defined by expressing the force (or moment) acting on the immersed 
surface associated with the dynamic response in terms of these 
parameters. Equation (4.10) gives: 

3 

6 

where when (i = 1,2,3) represents a force and when (i = 4,5,6) 

represents a moment. The parameters X_. and X_. denote the 
first and second time derivatives of the dimensionless dis- 
placements : 



j 1,2, ... 6 (4.12) 



\ . = - { 1 ;°Upa 4 M..X.+ poa 4 N..X.1, {. ’’ 

ij \ a M iJ J iJ jJ ’ i = 4 > 5 , 
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X . = X . /a 
J J 



and 



X — Xj / a , j — 1,2,3 



(4.13a) 




and 




j = 4,5,6 



(4.13b) 



5. NUMERICAL SOLUTION 



The primary objective at this juncture is to solve the 



integral equation (3.11) which has been established for the 



source strength function, f^(£,U»C). Once this function is 
obtained, u^ can be determined from (3.1) by evaluating the 



surface integral. Then, all other physical quantities such 
as pressure and resulting forces and moments may be determined. 

A numerical procedure can be devised by approximating the 
actual immersed surface of the structure by a contour composed 
of a large but finite number of facets. In the limit as the 
number of facets increases and the size decreases, the approx- 
imate contour approaches the actual contour. Thus, we may assume 
that in the limit as the number of facets approaches infinity, 
the numerical scheme converges. 

To begin, the immersed surface is subdivided as indicated 
in Figure 3 and an index is assigned to the nodal point (centroid 
point) of each of the small facets. It is required that (3.11) 
be satisfied at each of these nodal points. This is a relax- 
ation of the requirement implicit in (3.11) that the equation 
be satisfied at every point on the immersed surface. 
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FIGURE 3 NUMERICAL SCHEME 
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Recognizing that the source strength, f_.(£,n,4), is a regular 
well-behaved function, we may approximate (3.11) by 

-£ n u., y .,z.) + E i f „< x r y r z j ) ff*s j If (x i >y i’ z t ; x ’ y ' z) ds 



’ 2g n (x i’ y i’ Z i ) 



(5.1) 



where AS. denotes the area of the i-th facet and N denotes the 
J 

total number of facets. In shorthand notation, (5.1) may be written 



-f + f a.. = 2 g (5.2) 

n . n . xi n . 

x j J x 

where n = 1, 2, . . .7 corresponds to the six degrees of freedom and 

scattering, the repeated index indicates summation as usual and a 
is defined as 



ot . . 
iJ 





x,y,z) 



dS 



(5.3) 



Eq. (5.2) is a complex matrix equation which may be solved using a 

digital computer to obtain the source strength f once the elements 

i 

of the square matrix a., are calculated. 

ij 

As an approximation, the matrix a_ may be evaluated by using 
the value of 3G/3n at the nodal point to represent the mean-value 
over the elemental surface area of the facet and thereby further 
approximate (5.3) by 



“ij ■ 2? If (x i’ y l’ z i ; x ’ y ' z) 4S j 



(5.4) 



where 3G/3n may be evaluated by use of either (3.6) or (3.7) . The 
choice between the two forms depends on the value of r. When r is 
small, the infinite integral occurring in (3.6) is fairly rapidly 
convergent while many terms of the series indicated in (3.7) are 
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required. However, when r is large, the series occurring in (3.7) 
converged rather rapidly and, consequently, only a few terms are 
needed. 

The potential function u . may be obtained from (3.1) once the 
source strength f_. is determined at the nodal points. In a manner 
similar to that used to solve the integral equation (3.11), the 
source strength is taken outside the surface integral and (3.1) is 
wr i 1 1 en as 



U (x . 
n l 



z i } 



i 

4 TT 



N 



1=1 



f (x. 
n J 




x , y , z) dS 
(5.5) 



where as in (5.1) n = 1,2, ...7, and the integer N denotes the total 

number of nodal points on the immersed surface. In indicial notation 
(5.5) becomes 



u = f 3 . . 

n . n . ij 

i J 



( 5 . 6 ) 



where the square matrix g is defined as the integral over the 

panel of area AS. as 

J 



Ij = h fL. *>y. 



z) dS 



(5.7) 



Taking the value of G at the nodal point to represent the mean-value 

over the facet, (5.7) may be approximated for large R as 

= 7“ G ( x . , y . , z . ; x.,y.,z.) AS. (5.8) 

ij 4 tt i i i j j j 

where G is evaluated by use of either (3.3) or (3.4), depending on 

the value of r. When R is not large (5.7) must be used and the 

integration of 1/R term in G is discussed in Appendix B. 

The first step in the numerical solution to the problem is to 

ev aluate the matrices, a., and 3..* For example, the point ( x , y , z ) 
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on the immersed surface is denoted by i and the point where the 
source is located (5,n,£) is denoted by j. Thus, using the series 
form of the Green's function given in (3.4) or its derivative given 
in (3.7), there is no difficulty in numerical evaluation once the 
location of the point i and j are specified. However, the integral 
form of G and its derivatives as specified by (3.3) and (3.6) re- 
quires some special consideration. The infinite integral which 
occurs in both of these expressions poses two difficulties with 
respect to numerical evaluation; the integral has an infinite 
upper limit and the integrand is singular at y = y^ where y^ is the 
root of 

y tanh(y h) - v = 0 (5.9) 

o o 

The root is simply equal to "a" in view of (3.3f). However, these 
difficulties can be overcome by recognizing that the integrand is 
singular like l/(y-y^), subtracting the singularity out of the 
integral and carrying cut its integration analytically. The re- 
mainder of the integral is then carried out numerically and the 
upper limit is replaced by a suitable large number such that con- 
vergence is assured. 

For purposes of illustration let the infinite integral in 
(3.3) be denoted by I so that 

I - P.v/- (5.3. 0) 

J o ( y-y q ) 

wh er e 

P(y) (y-y ) 

o 

y tanh (yh) - 



Q ( y) 



v 



(5.11) 



and in the case of (3.3) the function P(p) is defined as 

(p+v)e ph cosh [ p (h+n ) ] cosh [ p (h+y ) ] J. (yr) 

P<l,) cosh(Bh) 



(5.12) 



The integral in (5.10) may be rearranged and written in the 
f o rm 



I = f 2 \ Q(y) ~ Q(y o ) dp + Q(p )P.V./* 

J O t N ° J 



2 p. 






dp 



0 (p- p Q ) 



(5.13) 



+ f C ° Q(y) d 

J (y- v Q ) 

where the principal value of the second integral can easily 
be shown to be zero and the other two integrals are proper 
and easily evaluated by numerical integration. Thus, I may 
be evaluated by numerical integration of 



i -/ (q(p.).-Q(u 0 ).) av+ f“ smih 

° (y - y o ) J 2p Q ( ^o } 



(5.14) 



where Q(y Q ) represents the limit of (5.11) as y-^y^ where 
y^ = a. The result of this limiting process results in the 

following expression: 

q( u ) = LLLcJ 

o' tanh(y h) fl-y h tanh (y h)"|+y h 

o l o o -J o 



(5.15) 



or, using (3.3f) and 



the fact that y =a, (5.15) becomes 

o 



n / _ a P ( a) 

P o v+(a 2 - v 2 ) h 



(5.16) 



For the case of deep water (h+°° ) we find v-^a so that 

Q ( y Q ) = P ( a ) 



(5.17) 
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The same method of integration may be used to evaluate the 
infinite integral in the derivatives of G as well as in G itself. 
The only difference is that P ( y ) is defined differently for the 
other cases but can easily be determined by a comparison of 
(5.10) and (5.11) with the particular infinite integral to be 
evalua t ed . 

An alternate method of evaluating the integral form of 
the Green’s function may also be used. This procedure as 
described in Appendix C involves converting the infinite upper 
limits and appears to be particularly useful when the wave 
period is large. The computer program actually uses this form 
for numerical calculation. 

There is one further significant difficulty in evaluating 

a . and g. .• When point i is distant from the source located 
ij ij 

at point j it is adequate to use approximates such as Eq . (5.4), 

and (5.8). However, when i = j or point i is near point j the 

singular term in the Green’s function which is of the form 1/R 

does not vary slowly over the panel of area AS^. Thus, when 

the node point designated by the index i is rather close to 

the j-th node point the 1/R term in G and derivatives of 1/R 

which occurs in a. must be integrated over AS. rather than 

ij J 

approximated. This integration is described in detail in 
Appendix B . 

Specification of the Subdivision Scheme 

The numerical scheme starts with specifying panel corner 
points on the immersed surface. Each corner node point is 
assigned an index and then a correspondence table is introduced 
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which specifies which four corner node points make-up a given 
panel. Thus, an index is assigned to each panel. 

Given the coordinates of the four corners of a given panel 
it is rather easy to compute its area by breaking it up into 
two triangular areas. Using the areas and centroids of the 
two triangular portions of the polygon the centroid of the 
polygon is computed. Finally, the unit normal vector for the 
panel is determined by taking the cross-product of vectors 
running across the two diagonals of the polygon panel. These 
calculations are described in detail in Appendix A. 

6 . DYNAMIC RESPONSE FOR A FLOATING BODY 

In this section the equations of motion for a floating 
body in waves are developed. These equations are then applied 
using hydrodynamic coefficients, which include the wave excita- 
tion forces and moments as well as added mass and damping co- 
efficients, in order to compute the dynamic response. 

The problem under consideration is represented schematically 
in Figure 1. In the following, the equations of motion are 
written with respect to the center of gravity. Thus, the origin 
of the body coordinates is assumed to be located at the center 
of gravity of the floating structure and d is defined as the 
dimensionless depth to that point. 

The small amplitude displacement of the body center of 
mass with respect to its mean position in the inertial reference 
frame is described by the three coordinates X^(t), X^(t) and 

X^(t) which are referred to as surge, heave, and sway, respec- 
tively. The small angular displacements of the body about the 
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x', y' and z' axes are denoted by 0 , 0 , 0 and are referred 

4 j o 

to as roll, yaw and pitch, respectively. 



The equations of 


motion linearized 


with 


respect to the 


small angular displacements of 


the body 


may 


now be written as 


f o 1 lows : 












T 

(t) = 


mX 1 (t) 








(6.1a) 


T 

F 2 ( t ) = 


mX 2 (t) 








(6.1b) 


T 

f 3 (t) = 


mX 3 ( t) 








(6.1c) 


T 

F 4 T (t) = 


1 , , 0, - 
xx 4 


I , , 0, - 

x y 5 


1**0, 
x z 6 




(6 .Id) 


F 5 T ( t ) = 


l , , 0, - 

y y 5 


l • » 0, - l , 

y z 6 y 


x' °4 




( 6 . le 


F 6 T (t) - 


i , , 0 , - 

z z 6 


I " 0. - I 

zx 4 z 






(6 .If) 


T 

where 


(t), F 2 T (t) 


and F 3 T (t) 


denote 


th e 


three components of 


the total 


external force acting 


on the 


body 


and F^ T ( t) , F 3 T ( t ) 


T 

and F^ (t) denote the three components 


of the total external 


mome n t . 


The symbol 


m denotes the body 


mass 


which equals the 


displaced 


mass . The 


moments of 


i ne r t i a 


are 


defined, typically 


as 


I , 

x y 


’ = I y ' z' 


f x T y T 

J m 


dm 


(6.2) 



where the integration is to be carried out over the complete 
mass of the body. For bodies having symmetry with respect to 
the (x f -y f ) and (y f -z T ) planes, all of the products of inertia 
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vanish. Although this type of symmetry is common to most ocean 
structures, there is no need to apply the limitation at this 
point since the inclusion of the product of inertia terms do 
not, in principle, complicate the development. 

At this juncture it is convenient to place the equations 
of motion in dimensionless form and replace the x, y, z sub- 
script system with a less cumbersome indicial system. We may 
define the following dimensionless parameters for this purpose: 



where, as previously given, a denotes the characteristic dimen- 
sion of the body, P denotes the fluid density and g denotes 




_3 



m 44 T x'x' /pa » m 55“ I y , y' /pa * m 66 -I z'z ,/pa 

m 45 =m 54 = - I x*y* /pi5 ’ m 46 =tn 6rt'i' ' ’ 5 

X 1 (t)-X 1 (t)/5, X 2 (t)-X 2 ( t >/a , X 3 (t)-X 3 (t)/a 

x 4( t > =0 4 ( t ) » X 5 (t)=0 5 (t), X 6 (t)=0 6 (t) 

f l T(t) = F l T( t )/pg a^ f 2 T (t)=F 2 T (t)/pgi 3 , f 3 T (t)=F 3 T (t)/pgi 3 
f 4 T ( t > =F 4 T ( t )/P8 s4 » f 5 T (t)= F 5 T (t)/pg5 4 , f 6 T (t)=F 6 T ( t )/pgi 4 




the gravitational constant. 



Using these definitions, (6.1) may be condensed to the 



f o rm 




(6.4) 



where ij takes on values 1,2, ...6, and the repeated index denotes 
summation as usual. 

Equation (6.4) represents the six equations of motion with 



coefficients as defined by (6.3). For free-floating bodies these 
coefficients represent the contributions from the surrounding 
fluid only and are generally considered to be composed of three 



dynamic forces and moments caused by the motion of the body, 
and (c) the hydrostatic forces and moments caused ^y the dis- 

k, 

placement of the body. For moored bodies the forces and moments 
associated with the mooring lines must be included and for the 
linear problem these contributions may be determined separately 
and superimposed. 

In accordance with this idea the three contributions to 

T 

the force (or moment) coefficient f (t) may be expressed as 
the linear combination: 




parts: (a) the wave excitation forces and moments, (b) the 




6 

c . ( t ) 4- [ c . . ( t ) + k..(t) + k!. (t)] (6.5) 

1 4 — ' iJ 13 

j=i 
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whe r e 



c . ( t) = force or moment coefficient associated with wave 

o / 

excitation, F^(t)/Pga' (i=l,2,3); F^(t)/Pgi (i=4,5,6) 

c..(t) = force or moment coefficient associated with the 

ij 

dynamic response of the body, see Eqs . (4.8-4.10). 

k..(t) = force or moment coefficient associated with linear 

ij 

or angular displacement which results from hydro- 
static pressures. 

k* . .(t) = force or moment coefficient associated with linear 

ij 

or angular displacement resulting from elastic 
constraints (mooring lines) . 

The force ormoment coefficient c^(t) may be expressed 



as 

i 5 

c.(t) = n ° R [IcJe 1 1 e~ i0t ], i = 1,2, ..6 (6.6) 

i * e 1 

where i = 1,2,3 refers to the three components of the force 

and i = 4,5,6 refers to the three components of the moment. 

The dimensionless coefficients c.(t) are defined according to 

T 

the definition of f. (t) as the force made dimensionless with 

1 

_ 3 4 

Pga (i = 1,2,3) or the moment made dimensionless with p ga 

(i = 4,5,6). The frequency of the wave excitation is denoted 

by a, and 6^ denotes the phase shift angle of the i-th component 

of excitation force or moment in relationship to the incident 

wave. (All phase shift angles are measured as lag positive 

and in relation to the time that the crest of the undisturbed 

incident wave is at the coordinate origin.) 
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The magnitude of the complex excitation force or moment 
coefficient occurring in (6.6) is defined, as in Section 4, 
a s 




Fj (max) 
-3 o 

pga n 

* 



i-1,2,3; | C i | 



F. . s 
i (max) 

-4 o 
pga n 

* 



i=4 ,5,6 



(6.7) 



where n° = H/2a denotes the dimensionless wave amplitude, H 
* 

being the wave height and a the characteristic body dimension. 

The amplitude of the excitation forces and moments are denoted 
by F-[_( max ) where i = 1,2,3 refers to the three components of 
force, and i = 4,5,6 refers to the moment. 

As the body responds to the wave excitation, dynamic 
pressures arise due to the motion which may be resolved into 
two components of force (or moment) , one component in phase and 
proportional to the acceleration of the body and a second in 
phase and proportional to the velocity. Equation (4.10) defines 
this dimensionless force (or moment) which is due to the motion 
of the body. The dimensionless parameter c_(t) is defined 
in (4.10) as the force (i=l,2,3) or moment (i=4,5,6) which 
is caused by the j-th component of motion. The dimensionless 
parameters M_ and N_ denote the added mass and damping coeffic- 
ients which are calculated by use of (4.6) and (4.7). 

The final contribution to the force (or moment) resulting 
from the surrounding fluid comes from the hydrostatic pressure. 

As the body is displaced from its equilibrium position, forces 
and moments arise which are proportional to the body displace- 



men t . 



The hydrostatic pressure increases with depth according 
to P = -P gy and, consequently, the i-th component of hydro- 
static force or moment resulting from this pressure variation 
and acting on the body is given by the integral 



F. 

1 



= Pga 3 



SL 



y g , d s 
s 1 



= 12 3* 



= Pga 



2,4 fL 



yg ± ds , 



i =4 , 5 , 6 



( 6 . 8 ) 



where y = y/a denotes the dimensionless y coordinate of a 

_ _ 2 

point on the immersed surface, and dS=dS/a , dS being an 
elemental surface area. The functions g ^ occurring in Eq . (6.8) 

are defined as follows: 



g-, = n 

1 x 



» = n 

2 y 



’3 = n z 



(6.9) 



g . = y f n - z 1 n , 

& 4 J z y’ 



g = z n -xn , 
s x z 



g = x n -y n 
y x 



where the dimensionless body coordinates are defined as 

x = x f /a, y 1 = y'/a, z 1 = z'/a, d = d/a 

and the unit normal vector on the immersed surface is defined 

as n ( x ’ , y 1 , z 1 ; ? , X ^ , X ^ ) = in^ + jn^ + kn^ . Eq. (6.9) is 

the same definition of g^ given in (2.21f), the only difference 

being that (2.21f) is written in terms of the fixed (x,y,z) 

coordinate system while (6.9) is expressed in terms of the body 



axes . 



The following linearized relationships exist between the 
coordinates of the fixed reference frame and the body coordinates 
for a given point on the immersed surface at x ? ,y' ,z' : 
x = x ' + X z T - X y f + X 



y = -d + y 1 + X^ x f - X.z ? + X, 

6 4. 



( 6 . 10 ) 



z = z T + X 4 y ? - X 5 x t + X 3 
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Equation (6.7) represents the hydrostatic force (or moment) 



acting on the immersed surface including the buoyant force 
which must, of course, just balance the weight. It is desired, 
however, to determine the i-th component of force associated 
with a j-th (j = 1,2, ...6) component of displacement of the 
body, j = 1,2,3 denoting linear displacement in the x f ,y f ,z f 
directions, respectively, and j = 4,5,6 denoting angular dis- 
placement about x' y y\z' axes, respectively. For the linearized 
problem this force (or moment) may be written, using Eq . (6.8), 



as 



H 



tj 3 f . _ _ 3 pi r r 

F_. X, = (1 or a) p ga X jf/yg. dS (6.11) 

J s 



1J 



9 X. j 
J 



(The factor 1.0 in brackets in (6.11) is applied when i = l,2,3 and 
a is appropriate when F^ denotes a moment (i.e.), i = 4,5,6.) 
Defining, further, the dimensionless parameter K as 



K ij f£ 8X 4 (y8 i 



) dS 



( 6 . 12 ) 



the dimensionless force (or moment) coefficient k. ,(t) may be 

ij 

written, in view of (6.11), as 

k^j(t) = -K^ X j ( t ) (no sum) (6.13) 

Equation (6.13) represents a form appropriate to Eq . (6.5) for 

the dimensionless force in the i-th direction caused by a dis- 
placement in the j-th degree of freedom. 



The following expressions are obtained for the hydrostatic 



force coefficients K. 

i J 



K 22 


* -ff s a y dS ’ V 3 * 


(6 .14a) 


K 24 


' K 42 -XC Z 'V S 


(6 . 14b) 


K 44 


■ /jC (y ' Z ' n z - z ' 2 ” y > ds 


(6 .14c) 


K 46 


■ K 64 ‘ffs*' z ' n y dS 


(6 . 14d) 


K 26 


■ K 62 - -ffs dS 


(6 . 14e ) 


K 66 


y ’n x - x ’ 2 n y> dS 


(6 . 14 f ) 



in which denotes the waterplane area in the equilibrium 

position. For floating bodies having symmetry with respect to 



the x ! -y T and y'-z f planes the coefficients, ^24 = ^42 = ^46 = ^64 



K 26 ' K 62 ‘ °' 

It is noted here that the expressions, Eqs . (6.14), for the 

hydrostatic force coefficients are in rather convenient forms 
for evaluation as a part of the same numerical scheme discussed 
in Chapter 5 which involves representation of the immersed sur- 
face by a large number of nodal points. At each nodal point it 
is necessary to specify the three coordinates of its location, 
the three components of the unit vector normal to the surface, 
and the elemental area of the facet. This information is, there- 
fore, available in convenient forms for use in the evaluation 
of the integrals in Eq . (6.14) in the computer program. 
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If the mooring lines can be approximated by linear elastic 



constraints then the force or moment coefficient may be expressed 
in terms of a spring constant K f _ and displacement in a form 

similar to Eq . (6.12) as 

k * ^ j ( t ) = -K* Xj(t) (no sum) (6.15) 

whe re 



i = 1,2,3; j = 1,2, ...6: force in i-direction/pga 



- 3 



k' (t) = 



i = 4,5,6; j = 1,2, ...6: moment in i-direction/pga 



- 4 



and K ! _ denotes a spring constant matrix which depends on the 
mooring line configuration, stiffness and tension. ^'ij 
de f ine d a s 





3F .1 

l 

3 X . 
JT 



_ 3 

/ pga 




/ pga 



i-1.2,3 

3 = 1 , 2,.. .6 

i=4,5,6 
j-1,2, . . .6 



( 6 .16a) 
(6.16b) 



in which F_^ = the force in the positive i-th direction caused by 
the mooring lines when i=i,2,3 and F'^ = the moment about the 
i-th body axis due to the mooring lines when i=4,5,6. The 
parameters X^ denote dimensionless displacements as defined by 
Eq. (6.3). That is, when i =1 , 2 , 3 > X^ denotes dimensionless linear 
displacemen ts 

X. = X. / a, i=l ,2,3 
ii 

and when i=4,5,6 angular displacements, 

X. = 0. , i =4 , 5 , 6 

i l ’ 
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The equations of motion for a free floating body may now 
be obtained by use of Eqs. (6.4), (6.5), (6.6), (4.10), and 

(6.13) as follows: 



— *" CJ * 

a (m. . + M . . ) X . + N . X. + 

~ J g 3 



(K . . +K ' . . ) X • 
ij ij J 




C . 
1 



e 



i 6 . 



-ia t , 
e j 



( 6 . 17 ) 



Furthermore, the body response in the j-th degree of freedom may 
be expressed in the form 

X.(t) = x! R [ e 1 e _iat ] (6.18) 

3 1 e 



which, when substituted into Eq . (6.17), yields the complex 

equations of motion, 



X C 

[- (m . ,+M. . ) -iN . . + (K . .+K ' . . ) / v]-r =1— l e l6 i 



(6 .19) 



ij iJ 



ij ij 1 3 



n. 



where the frequency parameter v= a 2 a/g. 

Equation (6.19) represents six equations corresponding to 
i=l , 2 , 3 . . . 6 . The repeated j index denotes, as usual, the summa- 
tion over the six degrees of freedom. The amplitude ratio 
X°/n* denotes the ratio of the amplitude of the motion to the 
amplitude of the incident wave and denotes the phase angle of 

the motion (displacement) in relation to the reference condition 
of the crest of the incident wave being located at the coordinate 
origin . 

It will be recalled that for bodies possessing x f -y 1 and 
y ! -z r plane symmetry, m^^O for i/j in Eq . (6.17). Also, M 

and N_ denote the thirty-six element added mass and damping 
t ens o r s . 
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It can be shown, however, that 



K . . = K . . , M . . = M . . , K . . = K . . 

ij Ji iJ Ji ij ]i 



( 6 . 20 ) 



so that only 21 different values exist in the general six degrees 
of freedom problem. 

To further simplify the notation, we may define the complex 
matrix , 



A . . 
iJ 



- (m. .+M. . ) + (K . . +K 1 . . ) /v-iN . . 
ij ij ij ij ij 



( 6 . 21 ) 



the complex response vector, 

„o 



= 



3 

r - e J 



and, in addition, the vector, 
i C . i 



B 



id . 

l 



Then, Eq. (6.1Sj)may be written very simply as 

A., y . = B. i,j = 1,2, ...6 
ij J i, J 



( 6.22 ) 



(6.23 ) 



Equation (6.24) may be solved by matrix inversion. Once the 
solution to Eq. (6.21) for is obtained, the problem is 
solved. The amplitude and phase angle of the response can 
be extracted from the definition of y^ given in Eq. (6.22;. 
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7 . HASKINDS RELATIONS AND ENERGY BALANCE 

Even though it may be supposed that the numerical 
solution proposed here will converge upon increasing the 
number of partitions, it is important to keep the partition 
size large (and number of partitions as small) as accuracy 
considerations will permit in order to reduce computer time 
and storage requirements. It is, for this reason, impor- 
tant to determine the effect of the partition size on accuracy 
so that practical limits may be established. 

One method of verifying the accuracy is to compare the 
numerical results with analytical results where closed form 
solutions exist. Although valid, this approach is limited 
to a few simple shapes; for more general shapes no such check, 
of course, exists. 

A second method of checking the validity of the numerical 
results involves the use of an energy balance as well as the 
use of the so-called Haskind’ s relations. Conservation of 
energy requires that a balance must exist between the energy 
required to oscillate the object, and the wave energy trans- 
mitted across some control volume surrounding the object but 
at a large radial distance. Using the asymptotic form of 
the Green’s function given in (3.4a) for large values of r, 
along with (3.1) the following relationship for the damping 
coefficient is so obtained. 
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N . . = 



a 2 - v 2 

1 77T 



ii 4 tt h ( a - v z ) + v 



cos (3 - 0) dS| 2 d0 



f Q 77 | J'J' f i (£ , n , c ) cosh [ a (h+n ) ] e lap l 

(7.1) 



where = /^^+r \ 2 and 3 = tan ^(c/5). This relationship 

expresses the damping coefficient in terms of the far field 
behavior of the solution. 

A relationship somewhat similar to (7.1) known as Haskind's 
relations, may be obtained for the wave force (or moment) 
coefficient. That is, the i-th component wave force (or 
moment) coefficient is related to the waves produced at 
infinity by the body oscillating in the i-th mode such that 



c i = cTo s h (ah') ffs f i u ' n, O cosh [ a(h+ n )] e ‘ x co s ( 3 -tt- ^ d g 

(7 . 2) 

where ip denotes the incidence angle as indicated in Figure 1. 
A form of this relationship between the wave force and the 
waves produced at infinity by the same body oscillating in 
otherwise still water was first derived by Haskind (14) and 
later reiterated and discussed by Newman (3). Equation (7.2) 
may be considered to represent a form of the Haskind 1 s 
relations as extended to the finite depth case. The details 
of this extension which involves integration by method of 
stationary phase are given by Rao and Garrison (15). 
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Equations ( 7 . 1 ) and (7.2) represent relationships for the 



damping and wave force (or moment) coefficient based on the 
behavior of the far field solution. A comparison of these 



pressure over the immersed surface, i.e., as obtained from 
(4.2) and (4.7), provides a convenient and valuable self- 
check on the accuracy of numerical results. These results 
are not limited to special configurations and may be applied 
to arbitrary shapes. Equation (7.2) is, however, limited to 
symmetry with respect to the x-y plane, a condition which is 
satisfied by most practical shapes. 

It may be noted, moreover, that for the special case of 
the vertical force on an axisymmetric body a very simple 
relationship between and N ^2 ma Y be obtained from Eqs. (7.1) 
and (7.2). That is, in view of the axisymmetry the integrand 
in Eq. (7.1) must be independent of 0. Accordingly, the 
integration with 0 may be accomplished by evaluating the 
integrands at any value of 0, say 0 = tt and multiplying by 2 tt . 
The right hand sides of Eqs. (7.1) and (7.2) then become 
similar and the following relationship may be derived between 
the wave force and damping coefficient in heave: 



results with N.. and C. obtained from an integration of the 

IX 1 




sinh (2 ah) 




(7 . 3) 



2ah + sinh (2ah) 
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For the case of infinite depth, Eq. (7.3) reduces to the 
special case presented by Newman (3) 




The numerical evaluation of Eqs. (7 .1) and (7 .2) can be 
easily carried out using the same numerical scheme outlined 
in Chapter 5. The numerical integrations over the immersed 
surfaces are converted to summations and evaluated using the 
source strength function obtained from solution of (5.2). 

8 . COMPUTER PROGRAM 

The theoretical development outlined in Chapters 1-7 has 
been coded for digital computer calculations and a description 
of that program is presented in this chapter. The general 
arrangement of the numerical procedure, a computational flow- 
chart and the input-output format are discussed. 

The basic requirements established for the program were 
considered to be as follows: 

(a) The program must work for floating bodies of 
completely arbitrary shape, but when either one 
or two vertical planes of symmetry exist the pro- 
gram should take advantage of this so as to reduce 
the input data and calculations. 

(b) Run time must be kept as small as possible but 
without losing accuracy. 
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(c) Computer core requirements must be kept to a 
reasonable value even for very complex shapes 
where a large number of panels are required. 

Flow Chart 

The flow chart for the computer program is given in Table 1. 
It should be noted that the total program is divided into two 
parts, a very small MAIN program which is used only for dimen- 
sioning, inputing parameter and calling the first subroutine 
DYNRES . This arrangement allows all of the subroutines which 
have variable dimensioning to be placed on files if desired 
and only the small MAIN program need be altered. 

The subrout ine OTECPR controls the flow of the calculations 
and calls essentially all of the other subroutines. It first 
calls the subroutine GDATA which in turn calls the subroutine 
CONFIG. This subroutine reads data cards describing the 
immersed surface and computes the unit normals, areas and 
centroid locations of each panel. GDATA then calls the sub- 
routine PLT which plots the panels on a calcomp plotter. The 
CPU time requirement to this point in the program is small, 
and therefore the program could be terminated at this point 
so that the plot output could be checked. Once it is deter- 
mined that the input data is correct a complete run could 
then be made . 

Next the subroutine WAVLIN is called by OTECPR . The 
purpose of this subroutine is to compute the value of "a" 
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TABLE I 



COMPUTER PROGRAM FLOW CHART - 0TECPR1 




OTTCPR 
Head . 
mass and 
snrinc , 
constant 
matrix 
cards 



r — \ 



Reads 
wave 
inci- 
dence 
angle . 



Dimensions 
Input parameters 




Plots the immersed surface on a 
Calcomp plotter. Input: ALF and 
BET, the azimuth and elevation 
angle of observer and R the dis- 
tance of the observer. 



Reads data cards describing 
the coordinates of the corners 
of the panels on the immersed 
surface. Computes the components 
of the unit normal vectors, 
areas of the panels and coordin- 
ates of the centroids of the 
panel s . 



WAVLIN 



Computes the value of a using 
the water depth and wave period 
in Eq. (3. 3f ) . 



GCOEF 



Computes parameters used in 

GREENS including the trial 

and error solution of Eq.(3.4b). 



0CINV2 


h 




| OCINV 


H 



Inverts ^--matrix 
(when NSS=2) 

Inverts <\- matrix 
(when NSS=1 ) 



out 

out 



of 



of 



core . 
core 



SCAT 



Computes <x>0 >cx x . c( r a z matrices 
and stores them on 'direct access 
files or in core depending on the 
value of the parameter, NS. 




Green’s function: series form. 
Green’s function: integral form. 

Computes the 1/R term in G. 

See Appendix B. 



VFT I Cnrinutf- 1 ^ p* . IIS id Vf j n bv T* f I . (3.21 f 1 
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Reads 
wave 
inci- 
dence 
anpcle . 



v 

0 



GCOEF 



• 1 1,-11 • • I i,l i 1 1 * ' . u •< (I 111 

GREENS including the trial 

and error solution of Eq.(3.4b). 



| 0CINV2 


h 


or 






| OCINV 


h 










SCAT 












| GREENS 


h 








or 


1 GREEN 







I GSING I 

~T , 

[ AVEVAL 



Inverts Qc-matrix 
(when NSS=2) 

Inverts (\-matrix 



Green's function 
Green's function 

Computes the 1/R 
See Appendix B. 



out of core. 

out of core 

matrices 
access 
on the 

series form, 
integral form. 

term in G. 



(when NSS=1) 

Computes Q( , ,(X x , 

and stores them on direct 
files or in core depending 
value of the parameter, NS. 



VELNOR 



Comoutes g^ as given by Eq.(2.21f) 




Out of core solution of Eq.(5.2) 

In core solution of Eq . ( 5 . 2 ) 

Computes the notential at the 
node points and prints the 
pressure and velocity dist . 



BAS COE 



HACK 



REG PON 



Using the nressure distribution 
BASCOE computes the excitation 
forces and moments, added mass 
and damning coefficients. 

Comoutes the value of the excita- 
tion forces and moments and added 
mass and damping coefficients by 
use of Haskind's relations and 
energy balance. 

Solves the equations of motion of 
the floating bodv,(i.e.), solves 
Eq. ( 6.24) . 



Loons through "N'VAVES" times. Each time through the loon a new 
wave incidence angle is read at the same neriod and denth. 
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by trial and error solution of the equation: 

v = a tanh (ah) 

The subroutine GCOEF is next called by OTECPR. The func- 
tion of this subroutine is to compute certain functions which 
are used repeatedly in the series form of the Green’s function. 
The first fifty values of p^ given by (3.4b) are computed and 
these values are used to compute the first fifty values of 
cos [y^(h+y) ] and s i n [ p ^ ( h+y ) ] for the various different values 
of y corresponding to the node points at the centroids of the 
panels on the caisson. The values of these functions are 
stored and are used by the subroutine GREENS for use in eval- 
uating the series form of the Green’s function. 

The subroutine SCAT is next called by OTECPR. This sub- 
routine evaluates the elements of the a and 3 matrices and 
stores them on FILE #10 and 14, respectively. SCAT calls the 
subroutine GREEN and G SING or GREENS to evaluate these matrices. 
GSING calls AVEVAL to evaluate the 1/R and 9(l/R)/3n term in 
the integral form of the Green’s function. 

The program has, in general, two options for solving the 
matrix equation (5.2) either by inversion of (a-I) through 
OCINV or 0CINV2 out of core, or by solution through COMAT in 
core. When the inversion out of core option is selected, 

SCAT calls OCINV or 0CINV2 to invert (a-I) stores (a-I) ^ in 

the file space originally occupied by (a-I). 
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In addition to the calculation of a and g matrices and the 



storing of these matrices out of core, SCAT also generates the 

derivatives of a (i.e., the matrices a , a and a and stores 

x y z 

them on Files 11, 12, and 13). 

At this juncture in the flow of calculations OTECPR reads 
the first of a series of data cards which specifies the incidence 
angle, ip , of the incident wave relative to the body coordinates. 
All of the computations made heretofore are dependent on the 
period and wave height but are independent of the incidence angle. 
Thus, computations may be made for a series of incidence angles 
at little additional expenditure of CPU time. The remainder of 
the program is, therefore, "looped" NWAVES times with different 
values for the incidence angle. 

The next subroutine called by OTECPR is VELNOR. This sub- 
routine computes the velocity induced normal to the caisson by 
the incident wave. More specifically, VELNOR computes the vector 
given by (2.21f). This vector represents the right hand side 
of Equation (5.2) so that the equation can be solved for the 
source strength. 

Depending on the particular option selected for solving 
(5.2) either SOURCE or COMAT is called next in OTECPR. The 
function of SOURCE is to multiply (a-I) ^ on g in order to 

evaluate the source strength through (5.2). In the case that 
(a-I) 1 has not been determined through OCINV or 0CINV2, the 

subroutine COMAT is called which solves (5.2) by elimination. 
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The subroutine BASPR is called next by OTECPR. The purpose 
of this subroutine is to compute the value of the potentials* 

(Uq + u y)» u 2 » u 2 9 u 3 5 u 4 5 U 5 an< ^ u 6 at the no< ^ e P°i nts ° n the 

immersed surface. Using these values of the potentials* BASPR 

computes the pressure distribution on the immersed surface. 

This subroutine also reads the FILES #11, 12* and 13 which store 

a , a and a and with this information computes the three 
x y z 

components of velocities in body coordinates at the node points 
on the immersed surface. It should be noted that the velocities 
are used only for a check on the solution and are not actually 
needed in the computation of the loads and response of the floating 
body. The computation of velocity can be surpressed by simply 
setting the parameter NBA = 1 in the MAIN program; when NBA = 2 
the velocities are computed. When NBA = 1 the FILES #11* 12 and 

13 are not used at all. 

The subroutine BASCOE is next called by OTECPR. The function 
of this subroutine is to re-compute the excitation forces and 
moments by use of the Haskind’ s relations (Eq. 7.2) and re-compute 
the damping coefficients, N^, by use of an energy balance (Eq. 7 .1). 
These results which are computed on the basis of the far-field 
potential are compared with the same results computed by use of 
the surface (near-field) pressure distribution as a check on the 
validity of the solution. 

Finally OTECPR calls the subroutine RESPON. This subroutine 
takes as inputs the excitation forces and moments, added mass and 
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damping coefficients, mass and spring constant information and 
then solves Equation (6.24) for the dynamic response of the body. 
Input Data 

The input data for the program describes the geometry of 
the immersed surface up to the mean waterline, specifies the 
water depth and wave conditions as well as the mass and mooring 
line spring constant parameters. Some of the data is input 
through data cards and some through constants in the DATA 
statements in the MAIN program. 

The 36 values of mass m and mass moments of inertia of the 
floating body, I-,-,, etc *> are input through six data cards, 
each card containing six values. 

The spring constants which account for the effect of the 
mooring lines are next input. Here again, there are, in general, 
36 values which are input through six data cards, each card 
containing six values. 

The immersed surface of the body is described through a 
set of data cards which give the three coordinates of the corners 
of the panels. These coordinates are measured with respect to 
a special coordinate system used for inputing the data only. 

These special ’’input axes" allow the user of the program to 
measure the coordinates of the panel corners in any convenient 
reference frame and then the subroutines GDATA and CONFIG shift 
the origin of this special input coordinate system to any desired 
location by use of the inputs XREF , YREF , ZREF and rotates it 
by the angle ANGREF. 
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Another set of data cards following the ones specifying 
the location of the corner points is used to specify which four 
corners constitute a given panel. Each data card in this set 
contains five integers, the first represents the panel number 
or index and the following four integers specify the indices 
associated with the corners. 

The last set of data cards specify the different values 
of the incidence angle of the wave. There are "NWAVES" of 
these cards . 

Coordinate Systems 

It is important to understand the coordinate systems used 
in order to input the program properly and, therefore, these 
are defined in the following. It is convenient to define 
three coordinate systems: 



(x" 


, y”, 


z") : 


Input axes 


(x' 


, y\ 


z’) : 


Body axes 


(x, 


y » z) 


: 


Inertial axes 



The input axes are used to input the location of the corners 
of the panels. This set of axes may be located at any convenient 
point and orientation with respect to the hull but the y"-axis 
is set vertical with positive y" measured upward. The second 
two coordinate systems also have their y-axes pointing upwards 
and the location of their origins and orientation is specified 
by the five parameters, XREF, YREF, ZREF, ANGREF and ECG. 
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The inertial axes is located such that the x ! -z ? plane 
represents the mean water level. It is located at the point 
x" = XREF, y" = YREF and z" = ZREF in the input reference frame 
such that the center of gravity of the body is on a vertical 
line passing through the origin. The body axes is considered to 
be aligned with the inertial axes when the body is in its equili- 
brium position except that it is shifted to the point (0, ECG, 0) 
in the (x,y,z) reference frame. When the center of gravity lies 
above the mean waterline ECG represents a positive number and 
when the center of gravity lies below the mean water level ECG 
will be negative. 

Figure 4 shows the three coordinate systems and their 
relationship one to another. The reason for the introduction 
of the input coordinate system is strictly for convenience in 
inputing data. In the case of certain configuration it may be 
convenient to measure the location of the corners of the panels 
in a coordinate system attached to the hull at a reference 
point which may differ from either the mean waterline or the 
center of gravity. 

Positive values of XREF, YREF, ZREF and ANGREF are shown in 
Figure 4. The wave incidence angle 9 \p 9 is also shown in the 
figure and is defined relative to the body and inertial axes 
as indicated. 

Corner Node Point Inputs - No Symmetry: 

It would be unusual that a floating structure would possess 
no symmetry whatsoever but the program has the capability of 
dealing with such general cases. When no symmetry is present 
the complete immersed surface of the hull must be covered with 
panels . 
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FIGURE 4 DEFINITION OF COORDINATE SYSTEMS 
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Corner Node Point Input - One Plane of Symmetry: 

When the immersed surface has one vertical plane of 
symmetry and the center of gravity lies in this plane, then only 
half of surface must be represented by panels. In such a case 
the input axes may be placed in any desired manner relative to 
the immersed surface but the parameters XREF, YREF, ZREF and 
ANGREF must be selected such that the plane of symmetry is coin- 
cident with the x f - y f and x - y planes and such that 0(x f ,y f ,z f ) 
is coincident with the center of gravity. The single symmetry 
case is depicted in Figure 5. 

Corner Node Point Inputs - Double Symmetry: 

If the hull possesses two planes of symmetry it is necessary 
to input data cards specifying the panel corners on one quarter 
of the hull only. Here again the input coordinate system i s 
used to input the coordinates of the corners of the panels and 
it may be positioned arbitrarily with respect to the hull. 

However, XREF, YREF, ZREF, ANGREF have to be input such that 
0(x,y,z) and 0(x f ,y f ,z f ) lie on the intersection of the two 
planes of symmetry and ECG must be such that 0(x f ,y f ,z f ) is 
coincident with the center of gravity. The result is that the 
(x,y,z) and (x f ,y T ,z f ) are placed relative to the floating body 
as shown in Figure 6. The portion of the immersed surface input 
falls in the positive x-z quadrant in the body and inertial 
axes as indicated in the figure. 
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FIGURE 5 SINGLE SYMMETRY 
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FIGURE 6 DOUBLE SYMMETRY 
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Although in the above two cases where the body possesses 
one or two planes of symmetry the input axis was defined such 
that XREF, YREF, ZREF, and ANGREF were non-zero, the usual case 
would be to input the data with input coordinate axes aligned 
with the body and inertial axes. This is generally convenient. 
However, it is not usually convenient to place the origin of 
the input axes at either the mean water level or center of 
gravity. It is generally more convenient to set the origin at 
the level of the bottom if the hull has a flat bottom so that 
YREF is normally non-zero. 

Input Fo rma t 

The inputs to the computer program are of three types: 

(1) specifying certain constants, (2) specifying the dimensions 
of the arrays, and (3) specifying the geometry of the immersed 
surface. These inputs and their format are discussed in the 
f o 1 lowing : 

Input Constants: 

The constants which are input by way of DATA statement 
cards are defined in terms of both the metric and English system 
in the following. It is noted, however, that either one system 
or the other may be selected but no mixing of systems is allowable 

PER Wave period and period of the motion of the body (sec 

H Mean water depth (ft. or meters) 

ABAR Reference dimension of the hull. Any value will 

do but normally the radius of the hull or the 
half-length is used (ft. or meters). 
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3 3 

RHO Density of the water (slugs/ft or kg/m ). 

2 

G Acceleration of gravity (32,174 ft/sec or 

9 . 8m/ s e c 2 ) . 



XREF 



Y REF 



ZREF 



x" coordinate of the origin of the body axes 
(measured in the input axes ) (f t . or meters). 
(See Figure 4). 



y" coordinate of the origin of the body axes 
(measured in the input ax es ) ( f t or meters) 
(See Figure 4 ) . 

z" coordinate of the origin of the body axes 
(measured in the input axes ) ( f t or meters) 
(See Figure 4). 



ANGREF . . . C lo ckwi s e rotation angle of the body axes 

(as viewed from ab o ve ) ( De gr e e s ) See Figure 4. 

ECG Elevation of the center of gravity of the body 

(y) measured in the inertial axes located with 
x-z plane on the mean waterline. 

ARMIN....Two different forms of the Green's function have 
been defined, an integral form as given by 
Equation (3.3) and a series form given by Equation 
(3.4). The series form requires much less CPU 
time to evaluate but the modified Bessel function 
of the second kind K (ar) becomes infinite as 
(ar)+0. Thus, when *(ar) is less than or equal to 
the value of ARMIN the integral form of the Green's 
function_is used. The parameter (ar) is given by 
ar = ( 27ra/ L ) ( r / a ) where L = wave length, a = ABAR 

= characteristic length scale of the hull and r 
denotes the horizontal distance between node points. 
A value of ARMIN = 0.01 is a reasonable choice. As 
ARMIN is increased the CPU time increases rapidly. 



RMIN The value of RMIN is used for two different tests 

of the distance between panels i and j in computing 
the values of a., and 8... The use of RMIN as well 
as ARMIN is besl:^ demonstrated by the following 
f low chart : 



7 0 



R < RMIN 




The value of RMIN is used to test R for purposes 
of selecting GREEN or GREENS. In subroutine GSING 
it is used to determine whether or not AVEVAL will 
be called to integrate the 1/R and derivative of 
1/R in the integral form of the Green’s function. 

The value of RMIN should be kept as small as 
possible consistent with adequate acc ur acy. Figure B1 
in Appendix B indicates that RMIN _>2 /AS* (where AS 
represents the area of a typical panel) represents 
a reasonable value fo RMIN. (ft or meters). 

NC Integer denoting the number of panel corner points 

on the portion of the hull which is input. This 
integer represents the actual number of data cards 
read by CONFIG which describe the location of the 
corner nodes on the hull. 

NB Integer denoting the total number of panels on 

the complete hull. The panels are assigned indices 
running 1 through NB . Figure 7 shows the numbering 
system when the hull has either single or double 
symmetry . 

NBA Integer having the value of either 1 or 2. If 

NBA = 1 the velocity components at the node points 
on the hull are no t computed and direct access 
files #11, 12 and 13 are not used. If NBA = 2 the 

velocity components on the immersed surface are 
computed. It should be noted that the velocities 
are not actually used in the calculations of the 
pressures, excitation forces and moments, added 
mass and damping coefficients or the resulting 
dynamic response. 
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DOUBLE SYMMETRY 



NUMBERING SYSTEM FOR THE CASE OF SINGLE AND DOUBLE SYMMETRY 
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NA Number of panels defined by the NC corner nodes 

on: (a) the complete hull for the no symmetry 

case, (b) half of the hull when the hull has one 
plane of symmetry, (c) one quarter of the hull 
when the hull has two planes of symmetry. In the 
case of (a), (b) and (c) NA will have the values 

NB, NB/2 and NB/4, respectively. 

NEQ The (a-I) matrix as indicated in Eq . (5.2) is a 

square matrix of size NB x NB in all cases except 
when NSS = 2, in which case it is NB/2 x NB/2. 

In all cases, however, the matrix is placed 
on direct access FILE #14 and is brought into core 
one piece at a time. The size of the block that 
is brought in is always NEQ (number of equations) 
by NB. Thus, it is most efficient to make NEQ as 
large as possible. However, the block size is 
limited (about 7200 bytes on an IBM 360) so NEQ 
may have to be small when NB is large so that the 
block size remains within the limits for the direct 
access device used. 

In general, NEQ may have any value equal to 
or greater than unity. The general rule is to make 
NEQ as large as possible consistent with the maximum 
block size allowable. 

NWAVES . . . In t eger denoting the number of wave directions 

which will be considered in the run. There must 
be NWAVES data cards giving the incidence angle 
(in degrees) of the particular cases to be calcu- 
lated. 

NSS Integer having either the value of 1 or 2. The 

normal value for NSS is 1 but under the special 

case where the wave direction is aligned with the 

x f -y f plane (0 or 180 degrees) and this plane 

represents a plane of symmetry, then NSS may be 

set to 2. The program will then take advantage 

of the fact that the source strength function on 

the positive z'-half of the hull is the same as 

at the corresponding points on the negative z f - 

half. Accordingly, the (a-I) matrix will be only 

(NB/2) x (NB/2) rather than NB x NB and, therefore, 

a savings in CPU time will be realized during inversion. 

NS Integer which is normally set to 1. However, 

under the special case where NSS is set to 2 the 
(a-I) matrix may be small enough to store in core. 

If this is the case, NS may be set equal to the 
integer NB/2. The (a-I) matrix will then be stored in 
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the ALS matrix and inverted through subroutine 
COMAT. It is always less time consuming to take 
the NS = NB/ 2 option when possible provided the 
storage requirements do not exceed that available 
since the in-core inversion is faster than the 
out-of-core inversion. 

Direct Access Files: 

The DEFINE FILE statements in the MAIN program specify the 
block size and number of records on the direct access files #10-14 
as follows : 



DEFINE FILE 10 (a,b,U,LK) 

DEFINE FILE 11 (a,b,U,LK) 

DEFINE FILE 12 (a,b,U,LK) 

DEFINE FILE 13 (a,b,U,LK) 

DEFINE FILE 14 (c,d,U,LK) 

where, for example, in #10, a = number of records and b = number 
of words in the record. The values of a and b must be input for 
a given set of data and conditions as follows: 

Number of records - 
a = NA 

c = NB/NEQ (rounded upwards to the nearest integer) 
Block size - (in number of words) 



b = 2xNB 
d = 2xNEQxNB 



The factor of 2 accounts for the 
fact that the words are complex. 



For the special case when NSS = 2 : 

c = (NB/2)/NEQ (rounded upwards to the nearest integer) 
In addition to the DEFINE FILE statements it is necessary 
to specify adequate and consistent disk space on JCL input. 



For efficient operation it is important that the five files 
#10-14 be placed on separate units when all are being used because 
these five files are written upon simultaneously in the same 
DO LOOP. 

When NBA = 1 files # 11, 12 and 13 are not used. In this 

special case it is possible to set the file space required for 
#11-13 to some nominal value. 

Summary of Program Options 

In view of the fact that the program has several options 
depending on the incidence angle and symmetry of the immersed 
surface, these are summarized in the following. The various 
options are completely independent of the values of NBA and NEQ. 

I. Hull has no planes of symmetry: 

NA = NB 
NS = 1 
NS S - 1 

The (a-I) matrix stored on file #14 will be NBxNB . The 
matrix will be inverted by use of an elimination scheme 
through subroutine OCINV. This option represents the 
most CPU time consuming method and would be used only 
if the hull had no planes of symmetry with respect to 
the body axes. The wave incidence angle may be arbitrary. 

II. Hull has either one or two planes of symmetry: 

NA = NB/2 (input data for +z f half of the immersed surface 
when hull has one plane of symmetry.) 

NA = NB/4 (input data for + x 1 and + z 1 quarter of the hull 
when the hull has two planes of symmetry). 



(a) NS 



1 



NSS = 1 

The program will take advantage of the symmetry in 
computing the matrices but the (a-I) matrix will 
be NB x NB and will be inverted in OCINV. The wave 
incidence angle may be arbitrary. 

(b) NS = 1 
NSS = 2 

The program will take advantage of the symmetry in 
computing the large matrices. The incidence angle 
must be set to either 0 or 180 degrees so that the 
wave is aligned with the x 1 -y 1 plane. The assumption 
is made that for this case the source strength on 
the -z T -half of the hull is the same as on the + z 1 
-half. Accordingly, the matrix stored on FILE #14 
will be NB/2 x NB/2 and will be inverted out of core 
by 0CINV2. 

(c) NS = NB/2 
NSS = 2 

The (a-I) matrix stored on FILE #14 will be of size 
NB/2 x NB/2. This matrix is eventually placed in 
core in the array ALS(NS,NS) which for this option 
has been dimensioned NS x NS, i.e., the size of the 
(a-I) matrix. The fact that NS is other than unity 
signals the computer to solve Eq . (5.2) in core by 

use of the subroutine COMAT. For this option 
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(NSS - 2) the incidence angle, ^ , must be 0 and/or 
180 degrees so that the wave direction is aligned 
with the plane of symmetry. 

Dimensioning 

The general principle to be followed in the dimension 
statements in the MAIN program is to make the dimensions exactly 
equal to the size of the array . In the case of two-dimensional 
arrays it is absolutely necessary that they be dimensioned 
exactly correctly. The dimensions cannot be either larger or 
smaller than the actual array. Even in the case of several of 
the one-dimensional arrays this is essential and, therefore, 
as a general principle the arrays should always be dimensioned 
to exactly the size of the array as given in the following with 
the exception of SH, CH, SINAMU and COSAMU arrays. 

The MAIN program contains all of the dimension statements 
which must be altered. It is unnecessary to make any changes 
in the subroutines because all arrays are variable dimensioned. 
The arrays in the MAIN program have the following dimensions: 
ALPHA (NEQ , NB) 

ALS (NS, NS) 

AV(NB), BV(NB), CV(NB), DV(NB), EV(NB) 

F ( NB , 7 ) , HH ( NB , 7 ) , U(NB,7) 

XB(NA), YB(NA), ZB(NA), XNB(NA), YNB(NA), ZNB(NA), DSB(NA) 
XC (NC) , YC (NC) , ZC (NC) 

XP(NC), YP(NC), ZP(NC) 

NODM (NA , 4) 

I C 0 M ( N S , 3) 

S H ( a ) , CH ( a ) , SINAMU(b), COSAMU(b), where 
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a >_ number of different levels of node points on the 
immersed surface. That is, "a" should be greater than 
the largest number printed under the column labeled LEVEL 
NO. on the output. A value of, say 20, would cover essen- 
tially all conceivable cases. 

b 5 0 x a (for example, b = 1000 if a = 20) 

(Note: In the case of SH, CH, SINAMU and C0SAMU the 

dimensions may be greater than or equal to the actual array 
size) . 

Input Data Cards 
Mass Matrix: 

The first set of data cards are used to input the mass 

matrix. There are six cards in this set; each card contains 

six numbers. The format for each card is 6F10.1 and the units 

2 

are slugs and slug-ft when the English system is used or kg 
2 

and kg-m when metric units are used. 

The first card contains the values of m 



11 



m 



12’ m l 3 ’ “‘14 ’ 



m. 



m n c and m- , the second card contains values of m nl m 00 , m 0 _, m_., 

1j Id 21 222324 

m _ _ , m 0/: , etc. where 

2j 2d 



m 



m 



m 



11 



12 



44 



m 
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m. r - I-,-, = moment of intertia about the z 1 axis (slues - ft 

6 6 z z . o N & 

or kg - m^ ) 



(Note 


, the x 1 , y ? 


and z f axes 


are centroidal axes 


by def ini tion) 
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Spring Constant Matrix : 

The second set of data cards are used to describe the 
elastic constraints of the body due to mooring lines. In the 
most general case, for displacement in each of the six degrees 
of freedom a force (i = 1,2,3) and moments (i = 4,5,6) can occur. 
Thus, in general, there can be 36 elastic coefficients to describe 
the effect of the mooring lines. 

The definition of the elements of the spring cons tan t 
matrix is as follows: 

3 F f 

rr I _ 1 



wh ere 



F . ' 

l 



X . 
J 



Examp le : 



Force in the i-direction (i=l,2,3)(lb or N) 
Moment about i-axis ( i =4 , 5 , 6 ) ( f t - lb or N-m) 



Linear displacement in the j - d i r e c t i on ( j = 1 , 2 , 3 ) 
(ft or me t e r s ) 

Angular displacement about j-axis (j— 4,5,6) 
(radians) 



The case of a simple spring restraint is indicated in the 



following sketch: 




It is possible to derive expressions for the values of K ! _ by 
starting with the expressions for the two components of force 
and mome n t : 



F 1 



1 



-k 1 (x 1 -x 1 ) 

o 



F ' 2 ‘ - k 2 < X 2- X 2 + l lV 

o 



-k 3 (X 2 - X 2 - <1 2 X 6 ) 

o 



F 1 & = -k 2 £ 1 (X 2 -X 2 + IX) 

o 



+ k 3 l 2 ( V X 2 ‘W 

o 



where the subscript "o" denotes the equilibrium position of the 

body. By differentiation, the spring constant coefficients 

K T . . are defined as: 
ij 
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where k^, k 2 , and k^ denote the spring constants of the three 

springs defined in the usual manner in terms of force/displacement 
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Hull Geometry: 

The next set of data cards to be read are the ones which 
describe the immersed surface of the hull. This set of cards 
gives the coordinates of the four corners of the panels and there 
will be NC of these cards. When there is no symmetry in the 
immersed surface, the complete hull must be described; when the 
hull has a single plane of symmetry (x ! -y ! plane), then half the 
surface must be described, and when the immersed surface has two 
planes of symmetry, one quarter of the hull must be described. 

When inputing half the immersed surface, the x ! -y ! plane must 
represent the plane of symmetry and when inputing one quarter of 
the hull, the x’-y* plane and y ! -z’ planes must represent planes 
of symmetry. In the case of single symmetry the +z’ part of the 
immersed surface is input, and in the case of double symmetry 
the +x f and +z ? quadrant is input. 

An example of a subdivided surface with two planes of symmetry 
is shown in Figure 8. On the quarter of the hull to be described 
there are 60 corner node points denoted and, therefore, for this 
example, NC = 60. It should be noted that all panels must have 
four sides; however, the length of the sides and orientation may 
be arbitrary, i.e., the panels need not be rectangular or square. 

The coordinates of the NC (60 in the example) panel corner 
points must be input by use of NC data cards punched in the fol- 
lowing format : 

Cols: 1-10: The corner node index (1-60 in the example) (Integer, 

right adjusted). 
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FIGURE 8 EXAMPLE BARGE CONFIGURATION 



Cols: 11-20: x M coordinate of the node (ft or ra). 

Cols: 21-30: y" coordinate of the node (ft or m) . 

Cols: 31-40: z" coordinate of the node (ft or m) . 

Correspondence Table: 

The next set of data cards will define the sequence of four 
corner nodes which make up a given panel. There are five integers 
on each of these data cards, the first represents the panel index 
and the next four integers correspond to the four corners of the 
given panel . 

There will be NA of these data cards where NA denotes the 
number of panels on the portion of the immersed surface described 
by data cards. NA will have the value: 

NA = NB (no symmetry) 

NA = NB/2 (single symmetry) 

NA = NB/4 (double symmetry) 

where NB denotes the total number of panels on the complete 
immersed surface. In the example described in Figure 8: 

NC = 60 

NB = 188 

NA = 47 

In this case the hull has double symmetry so NA = NB/4 = 188/4 = 47. 

Each of the NA cards will specify which of the four points 
compose a given panel and are punched according to the following 
format (all right adjusted integers): 
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Cols: 



1-10: Panel number (1 through 47 in the example) 



Cols: 11-20: First corner node index 

Cols: 21-30: Second corner node index 

Cols: 31-40: Third corner node index 

Cols: 41-50: Fourth corner node index 

The four corner nodes are specified in sequence as encountered 
in "walking" around the panel on the water side in such a way 
that the interior of the panel area remains on the right hand side. 
Another equivalent way of stating this is to number in sequence 
moving around the panel clockwise as viewed from the water side 
of the hull surface. The starting corner point is arbitrary. 

Wave Incidence Angle Data Cards: 

The final set of data cards are those specifying the different 
wave directions to be considered in the run. Since the (a-I) 
matrix needs to be constructed and inverted only once, it requires 
little additional CPU time to consider several values of incidence 
angle in addition to the initial one. It should be remembered, 
however, that when NSS = 2 the wave must be aligned with the x 1 
axis so the incidence angle can only be 0 and/or 180 degrees only. 

The format for the cards are: 

Cols: 1-10: Incidence angle, ip , in degrees. 

Computer Print-Out 

The computer print-out contains some of the input data as 
well as the computed results. 
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The coordinates of the corner nodes are printed as read-in 
except that they are shifted to the inertial axes (x,y,z) with 
origin at the mean water level on a vertical passing through the 
c.g. of the body. The units are the same as input. 

The correspondence table is also printed as read in for 
purposes of double checking. 

The centroids (node points) on the panels are computed as 
well as the components of the unit normal vectors and areas of 
the NA panels input. These results are printed under the heading, 
"COORDINATES OF CENTROIDS OF PANELS, COMPONENTS OF N AND PANEL 
AREAS. " 



The displacement of the immersed surface is computed by 
three different methods: 



xn dS 
x 



¥ -ffs 

V -ffs y " 

V ’ffs Z * 



dS 



dS 



The results for the volume computed in the three different manners 
represents a good quick check on the accuracy of the input data. 

If the input data is correct, the three results should be the same. 

In a similar fashion the total surface area as well as the 
area projected in both the x and z directions are computed as a 
quick check on the input data. 

Next the program goes through the NA panel node points and 
assigns an index to each different elevation. All of the node 
points at elevation -10.0, for example, may be assigned the index 
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1 and all at -5.0 would be assigned 2, etc. The panel node point 
index and its corresponding "level number" is, therefore, printed. 
Caution : largest number printed under "LEVEL NO." should no t 

exceed the dimension of the arrays SH and CH . The maximum value 
would typically be 3 to 10. A dimension of 20 to 30 should 
include all cases. 

The "WAVE NUMBER" is next printed. This is simply an index 
which indicates that a new wave direction is being considered. 

The number printed here will range from 1 through "NWAVES . " 

The next block of information in the print-out is the 
following : 

a = 27Ta/L = 2 tt . hull radius/wave length, 
v = (27r/T) 2 a/g 

h = water depth (ft or m) 

T = wave period (sec) 

ip - wave incidence angle (degrees) See Figure 1 
NB = total number of panels on total hull surface. 

The next block of data printed represents the pressure dis- 
tribution associated with the motion of the body in its six degrees 
of freedom: surge (1), heave (2), sway (3), roll (4), yaw (5), 

and pitch (6). The subscript (7) denotes wave interaction with 
the f ixed body . 

For each mode (six degrees of freedom 1-6 and 7 wave inter- 
action with the fixed hull) the dimensionless pressure amplitude 
coefficient PI, P2, etc. and pressure phase angle PHI, PH2, etc. 
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are printed. The definition of the dimensionless amplitude and 
phase angles are given on the print-out and will be given here 
also for completeness: 



PI = 


pressure (half) amplitude in height of water column /( ha 1 f ) 
amplitude of the motion in surge. 


P 2 = 


pressure (half) amplitude in height of water column/ (half) 
amplitude of the motion in heave. 


P 3 = 


pressure (half) amplitude in height of water column/ 
(half) amplitude of the motion in sway. 


P A = 


pressure (half) amplitude in height of water column/ 
(half) amplitude of the motion in roll times ABAR. 


P 5 = 


pressure (half) amplitude in height of water column/ 
(half) amplitude of the motion in yaw times ABAR. 


P6 = 


pressure (half) amplitude in height of water column/ 
(half) amplitude of the motion in pitch times ABAR. 



These definitions are the same as stated in Eqs. (2.20). 

The actual pressures at any given node point could be 
expressed as : 



p i = 


PI p gX x 0 a cos (PHI - ot) 


P 2 = 


P 2 p gX 2 ° a cos ( PH 2 - ot) 


U> 

II 


P 3 p gX 3 ° a cos (PH3 - ot) 


P A = 


PA p gX 4 ° a cos (PHA - ot) 


II 

PL| 


P 5 p gX 5 ° a cos ( PH5 - ot) 


P 6 = 


P 6 p gX 6 ° a cos (PH6 - ot) 


P 7 = 


P7 p g (H/2) cos (PH7 - at) 
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o 



denote the dimensionless half amplitudes of the motion 



where X. 

1 

as defined in Eq 



X °1 = X 1 °/a 
o 

X° 2 = X 2 /a 

X° 3 - * 3 °/a 

*°4 ' ®4° 

x ° 5 - e 5 » 

- °.° 



(2.1) where 



Since the motion of the body is defined by 

X. (t) = X°. cos ( <5 .-at) 

1 1 1 

it is clear the positive phase angles represent a lag of the 
pressure with respect to the displacement in any given degree of 
freedom . 

The next block of output represents the three components of 
fluid velocity (in body axes) associated with motion in the six 
degrees of freedom (1-6), and wave interaction with the fixed hull, 
(7 ) corresponding to each node on the complete immersed surface. 
(Note: if NBA = 1 this information is not computed and the 

output is suppressed) . 

The velocity component dimensionless amplitudes are defined, 
for example, as 

VX = velocity (half) amplitude in the x f -direction/ 
velocity (half) amplitude of the body. 

The phase angles are defined with respect to the displacement 
of the body (positive lag). The x’-component of velocity could 
be expressed, for example, as 
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surge 


(i) 


V =vx 

X 


X^ ao 


cos ( pxi - 


at) 


heave 


(2) 


V =vx 

X 


X° ao 


cos(PX2 - 


at) 








etc. 







The MASS MATRIX is next printed exactly as read-in through 

2 2 

data cards. The units are slugs and slugs-ft , or kg and kg-m . 

The SPRING CONSTANT MATRIX is printed next also as it was 
read in through the six data cards. Note that for a free-floating 
body all of the values in this array are zero. 

The next set of data represents the dimensionless wave load 
coefficients. The x,y,z subscripts represent body axes so FX, for 
example, represents the dimensionless excitation force in the 
x f direction. The definition of the coefficients and phase 
angles are given on the print-out. The actual wave force and 
moments can be expressed as a function of time as: 

F x = FX p ga 3 (H/ 2a) cos (PHASE FX-ot) 

F y = FY p ga 3 (H/2a) cos (PHASE FY-ot) 

F = FZ pga 3 (H/2i) cos (PHASE FZ-ot) 

M^ = MX pga^(H/2a) cos (PHASE MX-ot) 

M y = MY pga^(H/2a) cos (PHASE MY-ot) 

M z = MZ pga^(H/2a) cos (PHASE MZ-ot) 

Note: the wave crest passes the origin of the coordinate system 

at t = 0. Thus, a positive phase angle represents a lag with 
respect to the time the wave crest passes the origin. 

The HYDRODYNAMIC ADDED MASS COEFFICIENT MATRIX is next 
printed followed by the HYDRODYNAMIC DAMPING COEFFICIENT MATRIX. 
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These parameters are defined through Eqs. (4.9 - 4.13). However, 

to emphasize their relationship to the velocity and acceleration 

their definitions will be given here again. 

The force (or moment) in the i-direction (in body axes) due 

to motion in the j-direction is defined as F and is given in 

terms of the added mass and damping coefficients, M. . and N. as 

ij 1 3 

F.. = -(1.0 or a) ( p a 4 M.. X. + D o a 4 N.. X.} 



where 1.0 is used when represents a force (i = 1,2,3) and 

a is used when F_ denotes a moment (i = 4,5,6). In the above 
expression 



X. = acceleration in i-direction (Note, X. = X./a 
1 3 3 

for j = 1,2,3 and X^ = 0^ for j = 4,5,6). 

= velocity in j-direction 

The next block of output is labeled HASKINDS RELATIONS. 
Here the excitation forces, moments, phase angles are computed 
by use of: 



(1) the near field pressure distribution 

(2) the far field (radiation) potential by use of Haskind’ s 
relations 

and the damping coefficients are computed by 

(1) the near field pressure distribution 

(2) the energy radiated at infinity. 

The results computed by (1) and (2) are then compared and a 
percent difference computed. 

The use of the Haskind f s relations and energy balance 
represents a self-check on the accuracy of the results computed 
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by the pressure distribution. The results computed by use of the 
pressure distribution should be more accurate than those computed 
from the far field behavior of the potential but the two results 
should be fairly close. Large percentage differences indicate 
that there may be some problem and the solution may not be valid. 
However, if the values of the parameters being compared are very 
smalljlarge percentage differences should be expected and should 
cause no significant error in the final results. 

Finally the DYNAMIC RESPONSE results are printed. The 
response amplitude and phase angle in the six degrees of freedom 
are printed. The amplitudes are defined as: 



Surge: 


Half amplitude of the surge motion/wave (half) amplitude 


Heave : 


Half amplitude of the heave motion/wave (half) amplitude 


Sway : 


Half amplitude of the sway motion/wave (half) amplitude 


Roll : 


Half amplitude of the roll motion in radians x ABAR/wave 
(half) amplitude. 


Yaw : 


Half amplitude of the yaw motion in radians x ABAR/wave 
(half) amplitude. 


Pitch : 


Half amplitude of the pitch motion in radians x ABAR/wave 
(half) amplitude. 



The phase angles of the dynamic response are defined with 
respect to the crest of the incident wave in radians. The crest 
of the incident wave passes the origin at t=0 so a positive phase 
angle indicates that the displacement of the body lags the crest. 
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9. 



EXAMPLE COMPUTATIONS 



In this chapter a simple example calculation is carried out 
for a rectuangular barge (box) which is 90 x 90 meters in plan 
with a draft of 40 meters. Results are computed for a series of 
different wave periods. 

Hull Corner Node Points 

In Figure 9 is shown the subdivision of the hull; since the 
immersed surface possesses double symmetry it is necessary to 
describe only one quarter of it by use of data cards. 

The subdivision indicated in Figure 9 is rather coarse. 

There are 12 panels on the quarter-body and 48 panels all together 
on the total immersed surface. The quarter-body is described by 
use of the 19 corner node points so, 

NC = 19 
NA = 12 
NB = 48 

for the immersed surface described in Figure 9. 

The immersed surface was input with an input coordinate 
sy s tern placed relative to the body as indicated in Figure 9 
with the bottom represented by the y M = 0 plane. Since the draft 
was 40.0 m the inertial coordinate system was placed with origin 
at the position x" = z" = 0 and y M = 40.0. This is input to the 
computer program by use of a data statement specifying the 
parameters as follows: 

XREF = 0.0 
YREF = 40.0 
ZREF = 0.0 
ANGREF = 0.0 
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FIGURE 9 FLOATING BARGE 90m x 90m x 40m 



The computer "echo check" on the input data is shown in 
Table 2. The location of the inertial reference frame in the input 
reference frame is first printed in the table. The rotation angle 
is zero . 

TABLE 2 

IMITATION ANOLE= 0.0 



x -v : : f = 


0. 0 


Y - f = 


4 J. 0000 




O. 0 



r (v 


1 R : ) INA I lS 


JF l‘ANt : L CORN 


1 K 5 


u • 


X 


Y 


L 


1 


0.0 


-40.000 


0.0 


4 ? 


22.500 


-4 0. 000 


0. 0 


5 


45.000 


0.0 


0.0 


<* 


4 J . 0 0 0 


-20. 000 


o . q 


> 


45.000 


-4 0.000 


0.0 


6 


0.000 


-40.000 


22 .500 


f 


22.500 


— rC. 000 


22 . 500 


a 


*♦ 5 • 0 0 0 


0.0 


22.500 


9 


45. 000 


-20.000 


22 .500 


10 


4 5. 00 0 


-40. 000 


22.500 


li 


0.00 1 


0.0 


-* 5 . 0CC 


12 


22.500 


0.0 


5 • u o 0 


1 3 


45. 000 


0. 0 


v .> « o c c 


l't 


<♦5.0 0 ) 


-20.000 


4 5 . 0 0 C 


1 ^ 


<♦5.000 


-4 0.000 


-+ 5 .000 


16 


0.000 


-2C.000 


*♦ 5 . 0 C 0 


1 7 


22.500 


-20.000 


45.C0C 


id 


0. 10 ) 


-40.000 


4b . C0C 


1 ) 


22.50) 


-<t J. 000 


4 5 . 0 C C 
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14 
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6 
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1 
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5 


H 
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4 


9 


L l 
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1 7 


16 


iO 


12 
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14 


1 7 


Li 


16 




17 


19 


1 3 


12 


1 7 




14 


15 
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The coordinates of the corners of the panels are printed 
under the heading: COORDINATE OF PANEL CORNERS. A comparison 



of the numbers listed here can be seen to agree with the location 
of the corners (in the inertial reference frame) shown in Figure 9. 

Next, under the heading CORRESPONDENCE TABLE can be found 
the information which describes the given panels (numbered 1-12) 
in terms of the corner node point indices. Under NODE 1, NODE 2, 
etc. are listed the four corner node points which go to make up 
a given panel listed in clockwise order when observing from the 
"water side" of the panel. A comparison of these results with 
Figure 9 will confirm this. 

Table 3 gives the computer print-out of the location of the 
panel node points , the components of the outward unit normal 
vectors at the immersed surface (in body axes) and the areas of 
the panels. The node points are calculated by the computer, as 
described in Appendix A, and represent the centroids of the 
panels . 

In Table 4 is printed the volume and surface area of the 
immersed surface as computed on the basis of the information given 
in Table 3. The volume is computed by three different methods, 
i . e . , by 



n 




n AS. 
x . 1 



i 




AS . 

l 
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(z-proj ) 



n 

= L 

i=i 



AS 



The volume calculations in this way use all of the information 
in Table 3 and if the volumes computed by the different methods 
are equal to, say, 5 or 6 decimal places, it is quite likely 
that the input data is correct. This represents a simple quick 
verification of the input data but cannot replace the Calcomp 
plot for ease and reliability. 

The different indices assigned to the different levels of 
the panel node points are also given in Table 5. It may be seen 
by looking over the results printed in Table 3 that the node 
points fall at only three dif f erent levels: -10.0, -30.0 and 

40.0. As indicated in Table 5, and as can be verified by Figure 
9, panel nodes 1-4 are on the bottom and are located at the -40.0 
level (LEVEL NO 1), panel nodes 5,6,9,10 are located at the -10.0 
level (LEVEL NO 2), and panel nodes 7,8,11,12 are at the -30.0 
level (LEVEL NO 3) . 

The next block of computer print-out gives the WAVE NUMBER 
and lists the parameters, a = 27ra/L, v = a tanh(ah), where 
h = h/a, water depth, wave period, angle of incidence and the value 
of NB. 

Following this information, the computer print-out contains 
the pressure distribution and velocity distribution as previously 
described. After this an "echo check" on the MASS MATRIX and 
SPRING CONSTANT MATRIX is printed in the same units as input by 
way of the data cards. This information appears to be straight- 
forward and will not be given here. 



96 



TABLE 3 



uSDiNATrs or ccm r rcj t os nr panr ls , tc hx'n. ,\ rs o a»io p v;:l wpas 



X 


Y 


L 


NX 


NY 


N /. 
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44 9 . 4 99 
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- L 0. 000 
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0.0 


i .000 


4 4 / . 9 4 9 


3 3 • 1 5 0 
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4 3*00 0 


- 0 . coo 


0,0 


l .003 


4 4 ; * 3 'J 9 


U .2 30 
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♦ 4 . l) 0 0 
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0.0 


1 .cco 


4 4 V . 9 4 9 
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TABLE 4 

VUl.tX-PROJ.« 323999.125 

VOL. , Y -PR 0 J • ■= 323999.250 

V JL . » Z- p 4 )J. = 323999.063 



TOTAL SURFACE A«:;a= 22A99.961 

AREA, Y - (•' R G J . = -8099.964 

AREA, Z-F9 0 J . = 1799.997 



NUDE 



1 

2 

3 

4 

5 

6 

7 

8 
9 

LO 

11 

12 



TABLE 5 



level no 



1 

1 

1 

1 

2 

2 

3 

3 

2 

2 

3 

3 



0 .CO 
.0.00 
• 0 .CO 
,o.co 
LO. CO 

1 o . : o 

} 0 . vj 3 

30. CO 
LO . CO 
10. CO 
30. CO 
30 .00 



ELEV 
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The next block of information on the print-out is the dimen- 
sionless excitation force and moment coefficients. For the sample 
run this is presented in Table 6. It may be noted that the angle 
of incidence of the wave was zero so naturally Table 6 indicates 
FZ = MX = MY = 0.0. 

Table 7 contains the next three blocks of computer print- 
out, These three blocks represent the dynamic and hydrostatic 
force coefficient associated with the motion of the body. In 
making this example calculation, NSS was set to 1 even though it 
could have been set to 2 since the body is symmetrical about the 
x T -y ! plane and the angle of incidence was zero. Since NSS = 1 
all six degrees of freedom of the body were considered and, accor- 
dingly, all of the added mass and damping coefficients were com- 
puted. Had NSS been set to 2 only (M #1 , M._ and M . r ) and 

ll i2 16 

(N._, N . _ and N , . ) would have been computed. The value of the 
1 1 1 2 1 6 

other coefficients would be meaningless. 

It may be noted that the surge and sway added mass and 
damping coefficients, and and N^ and N^ ^ > are equal. 

This is, of course, due to the fact that the body is symmetrical 
in plan view. 

It is also interesting to note that the roll, yaw and pitch 

damping moments, i.e., N , N and N.., are quite small. This, 

4 4 j j ob 

in fact, appears to be typical of most cases, and, therefore, 
must be taken into consideration when interpreting that data. 

It may be recalled that the computed damping is wave damping 
only and this tends to be small in the angular modes. In such 
cases viscous damping may become important and it may be worthwhile 
to estimate the viscous damping. 
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TABLE 6 



**** WAVE LOAD CQEFF I Cl fcNTS AND PHASE ANGLES FOR THE LAISSCN 



FX, FY» FZ = AMPLITUDE GF FCRCS/R FC*G* A8AR**J (H/2 + A8AR ) 

MX, MV i MZ = AMPLITUOt OF M Ji'.ENT / R hQ *G * A BA R + * A ( H/ 2 * A B A R ) 

PHASE mNGLES = IN RADIANS MEASURED wl TF RESPECT TO TIP' 
WAVl CREST IS At ORIGIN (POSITIVE = LAG) 
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The next block of information printed in Table 8 is the 
Energy and Haskind's relations check on the accuracy of the 
solution. Here the excitation forces, moments and phase angles 
as well as the diagonal of the damping matrix is printed. The 
calculation of these parameters is carried out both through a 
pressure integration over the immersed surface and using the 
"far-field" potential through Haskind’s relations and the energy 
balance. The two sets of values and the percent difference 
between the two results are printed. 

It should be kept in mind that in this run FZ = MX = MY = 0 
so the percent differences in these cases is meaningless. Gen- 
erally, even for this very coarse grid the two sets of results 
agree rather well, at least when the values are not extremely 
small. For instance, FX and FY differ by only 5 to 6 percent. 

MZ, however, is rather small and the two results differ by about 
33 percent. This is, generally, the case; the moment tends to 
be rather sensitive to grid size and it requires a finer grid 
to get accurate moment results. 

The values of N ^, ^22 anc * differ by less than 6 percent 

when computed by the two different methods. However, the damping 
in pitch, yaw and roll are very small and tend to differ by a 
sizable amount. The difference is, however, not too important 
because such small values of the damping coefficient would not 
affect the response results except at resonance. 
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The dynamic response of the floating body is given in Table 
9. These results indicate that the surge (ha If -amp li t ud e ) is 
.5565(H/2) where H denotes the wave height (trough to crest). In 
this particular run the heave response shows a little dynamic 
amplification since the dimensionless response is slightly greater 
than 1.0. The heave half-amplitude is 1.2(H/2). The pitch half- 
amplitude is .0758 x H/(2a), in radians, where a = 45.0 m in this 
example problem. 



TABLE 9 

DYNAMIC RESPONSE 



Surge 


(X) 


** 


X/ETA = 


0.55650 


PHASE = 


1 . 51353 


Heave 


(Y) 


** 


Y/ETA = 


1 .19997 


PHASE = 


2 . 22306 


Sway 


(Z) 


** 


Z/ETA = 


0.00000 


PHASE = 


-1 . 11987 


Roll 


(X) 


** 


X/ETA = 


0.00001 


PHASE = 


-0.71514 


Yaw 


(Y) 


** 


Y/ETA = 


0.00000 


PHASE = 


0. 17472 


Pitch 


(Y) 


** 


Y/ETA = 


0.04758 


PHASE = 


-1 .61674 
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10. NUMERICAL RESULTS 



The purpose of this chapter is to present some typical 
numerical results for several simple geometric configurations 
and to make comparisons with experimental results. It appears 
that very little experimental data is in fact available for com- 
parison with the computed results but in cases where data is 
available comparisons are made. In addition, a comparison of 
the effect of subdivision size on typical results are presented. 
Numerical Results for a 90m x 90m x 40 m Free-Floating Box: 

Computer calculations have been carried out for the case 
of the floating box, 90 meters by 90 meters in plan by 40 meters 
draft, which was discussed in Chapter 10. The hydrodynamic para- 
meters as well as the dynamic response to wave motion are plotted 
as a function of wave period for two different grid sizes, 

NB = 48 and 108. 

The two different subdivisions of the floating box are 
shown in Figures 10 and 11 for the NB = 48 and 108 panel layouts. 
It may be noted that the figures show only one quarter of the box 
because the box has two planes of symmetry and, accordingly, in 
the two cases NA = 12 and 27. 

Model tests were carried out for the configuration in ques- 
tion at the Vassdrags Og Havnelaboratoriet (Rivers and Harbors 
Laboratory) in Trondheim, Norway. These results which were re- 
ported by Fa 1 1 i ns en and Michelsen (11) are compared with the 
numerical results presented here. 
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FIGURE 10: FLOATING BOX, 90M x 9QM x 40M. NB=48 PANELS 
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FIGURE 11: FLOATING BOX, 90M x 90M x 40M. NB=108 PANELS 
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The surge and heave excitation forces are shown in Figures 
12 and 13, respectively. The fact that there is very little 
difference between the two sets of results corresponding to 
NB = 48 and 108 indicates that the solution tends to converge 
very rapidly in increasing the number of panels (or decreasing 
the grid size), at least for the case of the excitation forces. 

It might be noted that slower convergence should be expected at 
small periods . 

The pitch excitation moment was computed also. However, 
as it turned out in this particular case, the moment with respect 
to the body axes (i.e., the centroidal axes was very nearly zero 
throughout the complete frequency range and the results were, 
therefore, masked by the "noise 11 caused by the numerical error.) 
For all practical purposes the pitch excitation moment was zero. 

The surge added mass coefficient is shown in Figure 14, and 
here again the results corresponding to NB = 48 and 108 indicate 
that the solution has converged. The experimental results pre- 
sented for comparison appear to agree rather well with the theory. 

The corresponding damping coefficient in surge is shown in 
Figure 15. These results are typical of most damping results. 

At high-frequency and low-frequency oscillation, the wave-making 
damping tends to vanish since no waves are generated at the two 
extremes. In the intermediate frequency range, however, the 
wave-making reaches a maximum and the damping coefficient like- 
wise shows a maximum. The experimental results shown on the 
figure appear to agree fairly well with the theory. 
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Fornax )/pga 3 (H/2a) 




PERIOD , sec. 

FIGURE 12 HORIZONTAL FORCE COEFFICIENT 
(90m x 90m x 40m, Infinite Depth) 



1 0 1 



Fy(max)/pga 3 (H/2a) 




PERIOD, sec. 

FIGURE 13 VERTICAL FORCE COEFFICIENT 
(90m x 90m x 40m, Infinite Depth) 
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FIGURE 14 SURGE ADDED MASS COEFFICIENTS 
(90m x 90m x 10m, Infinite Depth) 
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FIGURE 15 SURGE DAMPING COEFFICIENT 
(90m x 90m x 40m, Infinite Depth) 
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The added mass coefficient in heave, ^22’ s ^ own Figure 

16 along with experimental values. These results show very little 
frequency dependence, particularly when compared to the surge 
damping coefficient, N-^. This results from the fact that the 
surface primarily involved in heave force is the bottom and this 
surface is fairly distant from the free surface. Thus, free 
surface effects which are frequency dependent are not pronounced 
in the case of heave motion whereas they are much more pronounced 
in the case of surge where surfaces near the free surface are 
involved . 

The computed heave damping coefficient is shown in Figure 17 
along with experimental values. It is interesting to note that 
compared to the surge damping coefficient the heave damping co- 
efficient is rather small. This results from the fact that the 
bottom surface which acts as the primary wave-maker in heave is 
rather well removed from the free surface and is, therefore, not 
too effective in generating waves and resulting damping. 

It is also of interest to note that the experimental values 
of N ^ 2 fall well above the theoretical results. This difference 
is attributed to viscous damping and appears to be pronounced 
in a relative sense because the computed wave-making damping is 
small . 

The added moment of inertia coefficient shown in Figure 18 
is essentially independent of frequency (period). This appears 
to be typical of all results associated with caissons of the type 
in question. Moreover, the damping coefficient in pitch, 
which is not plotted, tends to be extremely small. These two 
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FIGURE 16 HEAVE ADDED MASS COEFFICIENT 
(90m x 90m x 40m, Infinite Depth) 
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FIGURE 17 HEAVE DAMPING COEFFICIENT 
(90m x 90m x 40m, Infinite Depth) 
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(90m x 90m x 40m, Infinite Depth) 
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features result from the fact that frequency dependence and 
damping is a result of wave production and very little wave gen- 
eration occurs as a result of pitching motion. This is due to the 
fact that the bottom surface is far removed from the free surface 
and the vertical sides do not move in the normal direction very 
much as a result of pitching motion. 

The added mass coefficient in yaw, M,.,., is shown in Figure 
19 along with experimental results. The agreement here appears 
to be generally good. 

Finally, the dynamic response of the box free floating in 

waves is shown in Figure 20. The pitch response,©,/ (H/2a) , is 

6 

extremely small as was expected. The natural frequency in heave 
occurs at about 16.5 seconds period and, therefore, considerable 
dynamic amplification occurs. In the case of the heave response 
several experimental values are shown which tend to agree well 
with the theory. 

Shallow*- Draft Barge 

The dynamic heave and pitch response has also been computed 
for comparison with experimental results for a shallow- dr af t 
barge. One quarter of the immersed surface of the barge is shown 
in Figure 21 as it was subdivided for numerical computations. It 
may be noted that since the draft was rather small the source 
strength might be expected to show considerable variation in 
the vertical direction and around the corner. However, in the 
horizontal direction the source strength should vary rather slowly. 
Accordingly, the panels were constructed with a rather large 
aspect ratio as indicated in Figure 21. 
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FIGURE 19 YAW ADDED MASS COEFFICIENT 
(90m x 90m x 40m, Infinite Depth) 
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FIGURE 20 SURGE, REAVE AND PITCH RESPONSE 
(90m x 90m x 40m, Infinite Depth) 
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FIG 21 ONE QUARTER OF IMMERSED SURFACE OF SHALLOW DRAFT BARGE SHOWING SUBDIVISIONS 



Calculations have been made to compare with the wave channel 
test results of Kim, Henry and Chou (17) for a shallow draft barge 
of 1.704 feet beam, by 2.558 feet length, by 0.16 feet draft. 

The location of the center of gravity and pitch gyradius is shown 
in Figure 22. 

The theoretical heave and pitch response is compared with 
the experimental values in Figure 22. The results show excellent 
agreement between the theory and experiment for the case of heave 
motion. This apparently results from the fact that the wave -making 
damping in heave for a shallow-draft barge is substantial and, 
therefore, the heave results are not too sensitive to viscous 
damping . 

The pitch results also show fairly good agreement. While 
there appears to be some scatter in the pitch results, it appears 
that near resonance the experimental values of the response fall 
below the theory. Here again this is attributed to viscous 
effects which tend to be significant compared to the rather small 
amount of wave-making damping. 

Disc Buoy 

As a further example, the dynamic heave and pitch response 
for a rather shallow-draft disc-buoy is shown in Figure 23. Un- 
like some of the previous results, the pitch response as well as 
the heave response shows excellent agreement between the theory 
and experiment. This is likely due to the fact that the corners 
on the disc buoy configuration are beveled rather than sharp. This 
probably reduces some of the separation and resulting viscous 
damping. Consequently, the experimental results are in good 
agreement with the calculations. 
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Floating Semi-submerged Sphere 



As an example showing the effect of finite depth on a floating 
body, numerical results are presented for a semi-submerged sphere 
floating in waves. Figures 24, 25 and 26 show the excitation 

forces, added mass and damping in heave, and added mass and 
damping in surge, respectively. In general, the effect of the 
finite depth tends to be more pronounced at low-frequency than at 
high-frequency. 

Using the results of Figures 24-26 the equations of motion 
for the free floating sphere give the heave and surge response 
presented in Figures 27 and 28. The amplitude ratios are shown 
in Figure 27 and the phase shift angles of the response measured 
with respect to the phase of the wave crest at the center of the 
body is shown in Figure 28. In general, there appears to be little 
effect on the finite depth for the sphere configuration. 

Floating Vertical Cylinder 

The added mass and damping coefficients for the circular 
caisson configuration are shown in Figure 29 as a function of 
the frequency parameter a ^ a/g wherein the depth ratio h = h/a, 
represents a parameter. The horizontal force coefficient is shown 
in Figure 30. These results for the horizontal force coefficient 
agree well with the theoretical results of Garrett (18). 

The dynamic response in surge, heave and pitch is presented 
in Figures 31, 32 and 33, respectively. In general, the results 

appear typical of deep water results and little effect of water 
depth is evident. The pitch response shown in Figure 33 shows 
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the large dynamic amplification at resonance which tends to be 
typical of essentially all floating caissons and, here again 
results from the rather small wave-making damping. This high 
resonant peak should not be expected in practice, however, due 
to the fact that viscous effects are undoubtedly considerably 
larger than the pitch wave-making damping used to compute the 
results shown in Figure 33. 
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FIGURE 24 WAVE EXCITATION FORCES ON A SEMI-SUBMERGED SPHERE 
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25 ADDED MASS AND DAMPING IN HEAVE FOR A SEMI- 
SUBMERGED SPHERE. 
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h/S=i.5 
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FIGURE 26 ADDED AND DAMPING IN SURGE FOR A SEMI-SUBMERGED SPHERE 



HEAVE AND SURGE AMPLITUDE RATIO 




FIGURE 27 HEAVE AND SURGE RESPONSE OF A SED.I -SUBMERGED SPHERE. 
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FIGURE 28 PHASE SHIFT ANGLES OF THE HEAVE AND SURGE RESPONSE OF 
SEMI -SUBMERGER SPHERE. 
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^IGU'E 29 ADDED HASS AND DAMPING COEFFICIENTS FOR A FLOATING 
VERTICAL CIRCULAR CYLINDER. 
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FIGURE 30 HORIZONTAL FORCE COEFFICIENT FOR A FLOATING VERTICAL 
CIRCULAR CYLINDER. 
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FIGURE 31 SURGE RESPONSE OF A FLOATING VERTICAL CIRCULAR CYLINDER. 
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FIGURE 32 HEAVE RESPONSE OF A FLOATING VERTICAL CIRCULAR CYLINDER. 




cr 2 a/g 

FIGURE 33 PITCH RESPONSE OF A FLOATING VERTICAL CIRCULAR CYLINDER 
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APPENDIX A: 



DETERMINATION OF THE CENTRO I DAL LOCATION, AREA AND UNIT NORMAL 
VECTOR FOR A PANEL 

The computer program is designed to accept the coordinates 
of the corners of the facets as inputs and, in addition, a set 
of cards indicating which corners make-up a given panel. A 

I 

panel is described by the four corner indexes read 
clockwise moving around the panel in such a way that the panel 
is always on the right. For example, panel 5, as indicated in 
Figure 1A, is specified by a data card 




FIGURE 1A 

of the following form: 

Panel No. Corner //I Corner // 2 Corner //3 Corner //A 
5 25 31 50 18 

In working with an individual panel the first corner (25) in 
the above example is designated as 1, the second corner, moving 
around the panel keeping the inside of the panel on the right , 
is labled 2, etc. The panel is re-labled as shown in Figure 2A 
with indices 1 -4 . 

The area of the panel is computed by first dividing it into 
two triangles by connecting corner 2 with the diagonal corner 4. 
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FIGURE 2 A 

The area of triangle 1-2-4 is then computed by use of a standard 
formula , 

Area (ii4) = ■ > /' 5 C s_ d lt )(s dz.<* )(S-W 

where 

S = 2^ C^'J- ^2.4- ^4( ) 

and 

d 1 2_ = ^J(x i_- x,) ~+ ( ~y z - y,) x 

dz.4- = - x t ) + ( % - yj 2- +• ( z 4 _ z 2-) 

d 41 = V(x, - x^) 2 - + (y,-y 4 ) 2 t Cz," 

The centroid of a triangle lies at the average of the co- 
ordinates describing its corners as, 

= -4f(*i + *>■ * 

y= = -j(v. + y* * y-o 

z <^ = i-c z ' + + *-**> 
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The area and centroid location of triangle 234 is computed 
in a similar manner. The total area for the panel is then simply 

A ~ **■ A n4 

and coordinates of the centroid of the panel is obtained from 



x c = 




A 12.4- 




y c ~ 


C Va 


A ,ZL^ 




A / (A i z.^ + A ^4' 


Zc = 




A JZ.-4- 


+ 


At.}*-')/ Az.3^- 



It is also necessary to determine 
the unit vector normal to the panel and 
the fluid. This is accomplished by defi 

■ P 

and r T t q. . which extend from point 1 to 
respectively, as indicated in Figure 3A. 



the x , y , z-components of 
directed outward into 
ning the vectors, 77 3 
3 and from point 2 to 4 , 
The corners are numbered 




FIGURE 3A 



clockwise when viewed from the outside (water side) of the caisson 
and, therefore, according to the right hand rule the cross product 

represents a vector having direction perpendicular to the plane 
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1 - 2 - 3 - 4 and pointing outward into the fluid. The unit 
normal vector is, then, obtained by normalyzing this vector to 
unit length, 



n = * % )/ I 'll, x % \ 

where jlft I = | and ft = L H* + 4 fl t , the 

n cj and n^. represent the direction cosines. 



components 



n, 
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APPENDIX B: 



INTEGRATION OF 1/R AND DERIVATIVES OF 1/R OVER A PANEL 

In order to evaluate the elements of the matrices Q( and [3 
by use of the definitions given in equations ( 5. 3^) and (5.7) 



integration is outlined. 

In order to carry out the integrations, a local coordinate 
system is defined in such a way that the panel lies in the x-y 
plane of the local coordinate system and the unit normal vector 
is parallel to the z axis. For definiteness the coordinate 
system is attached to the element in such a way that point 1 
is located at the origin and point 4 lies on the x-axis. 



the necessity arises to evaluate integrals of 1/R and 

of 1/R. In this appendix the method of 




y,°i 



PC*,q,2) 





FIGURE IB 



Consider first the integral 




ds 
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for the 



Following a procedure similar to Hess and Smith (16) 
velocity components, Faltinsen and Michel sen( 11 ) integrated 

equation (IB with respect to ^ to obtain 

-Jets y- -f- c*-*)*- +- ^ } 

s, 




The integrals in (2B) may be evaluated through numerical integ- 
ration in most cases. 



However, it may be noted that the integrand of the first 
integral is singular when z = O, J>= x and y - ^ <0. In 

this case the integral may be replaced by 



Jn[CX-S)^l- 






I 







(4B) 

d$ 
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The first integral in (4B) can be integrated analytically, the 
result of which is 



(x-s,) A [cx-$,)' u + - (/-sjA/'a-rjV z'-J 



(5B) 






£ z[ xJL'J - -ian’^Adb) )] + Z [(/- Zj - (*- S,)] 



The second integral in (4B) has no singularity in the integrand 
and consequently presents no dif f iculaties in numerical integ- 
ration. Difficulties with singularities in the integrand of 
the other integrals in equation (2B) are handled in a similar 
manner . 

In the evaluation of as defined by equation (4.5) or 

yC.. , Y u .- and Y .. as given by equations (4.11) it is 
lj "u *-j 

necessary to evaluate integrals of the form: 



AS 



(x-S ) 

Rj 






(6B) 






(7B) 



^ * fL ** 



cl S> 



( SB ) 



where 



R - V 7x^Ty~ 1 'c~y-Yy L + 



Equations (6B-8B) have been integrated by Hess and Smith' 
the results of which are summarized as: 
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(9B) 




A, -t-nx.- d rt. \ 

d/z. ^/i z + d iz ) 



?j ~ ~h Jin f h th Z-d+i 'j 
d Li / 



?4 ~ 7 3 A, X {La t. dh—z gLu£ ^ _i_ h ~ /jdL± th± <^±r_ i 

d 3 + ( ft -3 <-/L+ -hd 3 + J d+, \ /l ++A, * c 4 / / 




<=U 





■f /li - clti \ 


+ dz- 5 3 (L ( 


(• /*, 


+ Si t +d u J 


d, 5 ^ ( 



yi l f- /1 3 — c/o \ 

ft-z. + ''Z-s J 



(10B) 



i~-3~ 1±. JL ( di ± sl+ jz 4 j 5- ) + d LmMj Li f {hlAin 4±±. 

<=U M l * *+ -d» J d 4 , 



d) - j'a.n l . CHia e, - h L 
l Z A, 



_ /a n ( JUll^j- th 



+ 



j'an th. ' j _ iari 1 ^ 



2 /Z- 2 . 



i(xri‘ ( (dn ^-3 ~ th 
z a 3 



-j- -fan. Ldiidji — th j _ -far ? ^ (di±J^± L± ^ 



(HB) 



4 - Lari' f Jd±L^Jh) _ -fan 1 f m *‘ — L-Z—lh 

V 2 ^ J \ Z n, 



where 

d,i = /(?,- r,)V adw 1 

dii = V(S 3 - Oj- ^ l ) >_ 

^ = y><-- gy^-t- cu- ?,)*• <12B) 

d *. - yes. - 0 % ci.-%y 
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in which 



m,2. = 



h. 






7, 






^z 3 



la -T_.?_L 

§3-^ 



/77 3 ^ 



7 * - ?3 



/?? 



-4/ 



7, - Lfc 

5 ; - 



and 

/Ife = V^-5fe)^ *- (y - ? k ) i -(- > fc=/,z,J,4- 



( 13B ) 



(14B) 



e fe = Z 1 f Cx-t^ , M^/,2.,3^ (15B) 

V = (V-IOCX-S*.) > k-1,^3,4- (16B) 



In actually evaluating these expressions and cause 

no trouble. They become infinite on the edges of the quadrilateral, 
but in practice they are never evaluated there. The component 
requires special handling in certain cases. As 3Z o , 

O if the point P is approaching a point in the plane out- 
side the boundaries of the quadrilateral. If P approaches a 
point within the quadrilateral -*■ Z IT ( sgn z) as z -*■ o 
These facts may be verified from equation (11B). In the course 
of this method of flow calculation it is required to evaluate 
<P at points in the plane of the quadrilateral elements. In 
particular, the centroid of each element is in the plane of that 
element and within the quadrilateral. At such a point <p^ should 
equal Z'iT , since the case of interest is that for which z -* o 
through positive values, rather than through negative values, 
i.e., the field point approaches the caisson surface from the ex- 
terior flow field rather than from the interior of the caisson. 
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It may also be required to evaluate induced velocity components 
at points in the plane of a quadrilateral outside the bound- 
aries of the quadrilateral, for example, at the null points 
(centroids) of other elements if the body surface has a flat 
region. Points in the plane of a quadrilateral element should 
have z = 0, but, because of round-off error, they may have have 
small values of z with either sign. Thus, for points inside 
the quadrilateral, equation (11B) may give ZTf with either sign. 
To avoid this error, the absolute value of z is tested before 
velocities are computed, and if it is less than some small pre- 
scribed number, which is nevertheless large compared to the 
expected round-off error, it is set equal to plus zero and each 
inverse tangent of equation (11B) is set equal to Tr/2 with 
the sign of the numerator of its argument. This gives cp z — O 
for points outside the quadrilateral and cp> z - zrr for points 
inside the quadrilateral. Another situation that may cause 
trouble in the computing machine is when the slope of a side 
of the quadrilateral is infinite, i.e., when a side is parallel 
to the y-axis. It is evident from equation (11B) that in this 
case the two inverse tangents corresponding to that side cancel 
each other. To avoid difficulties each of the quantities 
( - % ) , ( J 3 - r L ) , ( y*. - 3Fj ) , and ( J ( - Xf.) are tested 

to determine whether they are zero, and if any one of them is 
zero, the two inverse tangents corresponding to that side are 
set equal to zero. Finally, it should be mentioned that the 
inverse tangents in equation (11B) are evaluated in the normal 
range -7T/2 to -t-Tr/2. It is tempting to combine some of the 
inverse tangents in this equation using the tangent addition 
law, but if this is done, great care must be exercised with 
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regard to the range in which the resulting inverse tangents 
should be evaluated. 

When the distance R is large it is possible to approximate 
integrals in Eq.(lB) and (6B-8B) by simply multiplying the integrand 
by the area AS. Thus, in order to get some indication of how large 
R must be in order for this approximation to be valid Figure (IB). 
The results shown here indicate that in the case of either the 
induced velocity or the potential contributed by the panels of area 
AS it is necessary that R2 2 / { a ~ s ' in order for the approximate 
expression to become valid. 
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FIGURE LB (a) COMPARISON OF EXACT AND APPROXIMATE ESPRESSIONS FOR INTEGRAL OF 1/R 




K-A-h as=ab=i 
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FIGURE IB (b) COMPARISON OF EXACT AND APPROXIMATE EXPRESSIONS FOR INDUCED VELOCITY 




APPENDIX C: ALTERNATE FORM OF THE GREEN'S FUNCTION CONVENIENT FOR 

NUMERICAL EVALUATION 

Except when the point x,y,z lies close to the source 

located at ^ , J the series form of the Green's function given by 
Eq.(3.4 ) is used for computer numerical evaluation. The series 
appears to converge fairly rapidly and, therefore, the numerical 
evaluation is quite efficient. However, when the point x,y,z lies 
close to the series form cannot be used because of the 

singular nature of the Bessel functions ">^,(d/i)and h^(OA.)as /t— »- O . 
In such cases, i.e., when ( Q/t) becomes less than some arbitrarily 
small number, the form of the Green's function specified by Eq.(3.3 ) 
must be used. 

However, Eq . ( 3 . 3 ) may not be evaluated numerically in a direct, 
straightforward manner because the 1/R term is singular as fc-*- O 
and the integrand of the infinite integral is singular at r. = a. 

The evaluation of the 1/R term through numerical integration was 
discussed in Appendix B. In this appendix an efficient method for 
the numerical evaluation of the infinite integral is discussed. 

The integral form of the Green's function is given by 

G(x,y,z;g,'} l l) = -jC +. -L 

(1-C) 

+ 2. Pv. COs hf-H(l+h)l ( M/t ^ djt 

J 0 -^< s/nJ? mIi - V cosh ,uh 

+ L ZJL C 2LD- cnsfiLa CLthll cosh Co. QdjA ll 3 * (ay-) 

a z h - y z h -r v 
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where 



R. — y (x^Ry- + ( </- > y + (z- ~s) z 
r 1 =~Y?x^s cy+zh-hj?' + (z- jy ' 

yz_ = yry^iy - ^ cz - xy 

y — cr z g _ q -{-arch (ah) 

The integrand of the infinite integral can be simplified by use of 
the relationships, 

swk(x) ~ cosh (X) ^ <2* 

for x larger than about 4.0. Thus, the interval may be broken into 
two parts. Using, for brevity 

/ \ _ Cm + \j 1 (£^L cosh RlCJiAIL cosh (y+h )] -Jf-n* ) (1 , c ^ 

' ' £'nh su/i — -y’ cosh /■ ih 



Eq.(l-C) may be written as 

r°° r u > 

2 RV. / ( ) dcL - Z RV. I / ) d y u 

do V ' 

do ( 3_c ) 

where U, 2. ^.o/h) . The term (m+Y)I(m-V) may be written 

^ / 

so that the infinite integral becomes 

XOO O U ! 

( j y = 2 pn ( 

Oo 

i pyj a-.(MX) dy ' 

J u, 

CO 

f /-?y f jrmsR r u 

J u M- v 



This may be further rearranged by use of the relationship given 
for example, by Gradsheteyn and Ryzhik (21). 
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_J_ — f e^^ ?? 32 (M A.) du 

P" 1 



(5-0 



where R" represents an image source with respect to the free sur- 
face 



R 



// 



— (. l j+ ' o) 



(6-C) 



The infinite integral may now be written, 



Z P.V.[ ( ) dyU = 



+ 2 J 



oo 

'' u ‘ e ^(*+}) j- d(M/L) d/4. -t PV.f 

J. yu-y 



(.v+V 



(7-0 



J C W d/i 

~U 

"O 

The remaining infinite integral can now be expressed in terms of 
a finite integral by use of the relationship (See Abramowitz and 
Stegun (22)) . 



dl(M^) = / f 



TT~ 



lm/lcos O , 

d& 



( 8-0 



Using Eq.(8-C), the last term may, therefore, be written as 

oo tj - 



,oo 



p. v -£±L e 

ja ~y> 



-j- 



u, 



u, L 



9-C) 



Now, let t = |t- y> and dt = dj 4 in Eq.(9-C): 

J 0 du - AUi 



u, 



TTs-bO 

y- fp -J j J r— ^ cte/ai 10-0 

° t-U-V 



or, simplifying further by use of the substitution S - t/(u^-v>) 
oo 

RV.J t Jd- v e J~o(MA) d/( 

cose) (u- »)s e vc ^Scho)#. 






r°° r 71 ' , 
= [ /_e v 

rr J J s' 

5=io 



d s de- 
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The exponential integral as defined by Abramowitz and Stegun is 

„co 






°o 



dt = 

j/ L 



e ^ 



( 12-C ) 



t 






so, the above result may be written as 

'' e ° 2» /=> M(LJ +V 



R V. 



'u, 



A - y 



J~ a (m/i.) cl 'a 



( 13-C ) 



, r (mi -*■ i a. cos q) 

2-R I (0 ' E'l-Cy+p ec Acos&) d& 



With this, the complete expression for G may be written 

At —h 

Q - -L- ,/./,* pi/ / ( - u+ rrshfM (^h)) Cosh[M(y-di)] U RR cUc 

2 IE' 12” J a Cosh/ih (M Unhph - p) ( 14-C ) 



Li ' 



ECmr.) dp 



o 




\>( j case) 



E'l-Cy-i-ii-i Acos&Xn-pj) de 



-h i ±KCaEi >21 cpsJiEl JlAll coshJAELLhil JniQdO 

Cdh - p L h t- R 



This form of the Green's function involves no infinite integrals 
and, therefore, can be evaluated numerically fairly efficiently. 



The integrand of the first integral is, however, singular 
when ja = a so a special procedure must be followed in evaluation 
of that integral. Let 



f(M) = 2 (Mi y) e M>1 Cash [a ( 3 + h)'} cosh [a ( yeh) 7 Jo (MnJ) 

Cosh (m h ) 



( 15-C ) 



The integral in question can be written 

2d 

r u> r ZQ r f 

RV\ C/A _ / -FCmI - U* L Ju + f(Q) —7 

J 0 M -fan h m h - y> J M (ann/Ah - y JM-fa 



2m 



m fanh mH - y 



A, 



1 

2d 



{jdii ddd 

M ~hxnh uh - R 
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The first and third integrals in Eq.(16-C) are not singular and 
will cause no difficulty in numerical integration provided the 
integrand is not evaluated exactly at the midpoint where yu = a. 

The integration procedure used, therefore, is the "3/8 rule" which 
evaluates the integrand at end points and two intermediate points. 
Using an odd number of subdivisions then the integrand is never 
evaluated at ^ = a. 



The second integral in Eq.(16-C) is singular and a special 
procedure must be followed. The integral may first be written as 



.24. 



d 



M 



* U 0 



where 



M ianh M-h ~ y 

Uo = ah 



J. 



U_ 



U i-anh a - Uo hnnh U 0 



( 17-C ) 



which indicates clearly the singular nature of the integrand at 

U-- U 0 . The interval may be broken into three parts as 

U 0 -£ 



f 



£Uo 



jDlLL 



c/u 



(18-C) 



J U ianh u - U c tanh U a ~ J a tank u - u e tanh u a 
r u 0 +a 

f du dJA 

t~ / 'll tanh u- Uo ianh iU ^ ~u tanh u - ~u 0 iduaFu 0 

U-o U'i-6 

and the first and last integral integrated in a straightforward 
manner by numerical integration. The second integral is evaluated 
by first expanding the denominator in a Taylor series about 
and then dividing 1.0 by this series. The result of this is 

u o) ~ Qh + ~d(~ah ~~ ^)( u ~ u °)"]da 



r da 

U tanh U- ll 0 tanh U 0 



11,-6 



Uo+<6 

= ffod 



U -£ 

O c. 
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— G 



a, a 



^ + ~L (a.± - Ik ) A + ...] d\ 
aj- a, l a f q, / J 



( 19-C ) 



= - % (2£) +- QCdE*) 

^ i 



where 



a l 



a 2 



i-an h (ah) + Q.h sech z (ah) 

sech\ah) ( !~ ah +anh (ah / 

= _L sech L (ah)[zah(\- Sech Z (*G )-- 3 hanh(ah)-ahsech z (ah)J 



Thus, the complete integral, Eq.(17-C) may be written 



ZA. 



du 



Uo~£ 

r 



m ianh/Ah - V 



du. 



-t- 
Lift 6. 



U f~anh U — fanh U 0 

(20-C) 

_ 2<£ sech z fah) jink gAJ 

u tanhu-u „ W,i/„ (Janh(ah) t ah scch\ah))^ 



ZUo 



dm 



where £ denotes some small value between 0 and Q.h . In the numer- 
ical evaluation, £ has been set to 0.10. 

The final form of the Green's function in the form appropriate 
for numerical evaluation may now be written as 



z<x 



U, 
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-f(a) f dff 



su fanh yu - ah -fanh (ah) 



zah 

u, + f(a f- —* 



c/u 



M fanh M- ah -fan hah) 
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+ zvfe 9i ^ +Uccs ^ E,L-( V +j + iAc°seXu,-v>] de 
Jq 

• ZJT_ (_£^Z Hl2. cosh faf 3 i-h )J Cosh£o-( h )j J^XOja} 

a l h - v-/? +. v 

where 

r. . _ 2 + y) cosh L m (1 + h)l Cosh f/U C y + h ) ] JZh' x ) 

tw - cos TTUT — 

The derivatives of G(x, y , z;j,^ , J ) follow from Eq.(21-C)as: 



M= - j- c*-s) - - ^3 
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- 4-.., (*-$) 
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where 

fw = 
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and finally, 
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